{'fraction_non_alphanumeric': 0.03988434620477207, 'fraction_numerical': 0.012673278028481097, 'mean_word_length': 4.424802110817942, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hybrid Particle-Continuum Simulations Coupling Brownian Dynamics and Local Dynamic Density Functional Theory

We present a multiscale hybrid particle-field scheme for the simulation of relaxation and diffusion behavior of soft condensed matter systems. It combines particle-based Brownian dynamics and field-based local dynamics in an adaptive sense such that particles can switch their level of resolution on the fly. The switching of resolution is controlled by a tuning function which can be chosen at will according to the geometry of the system. As an application, the hybrid scheme is used to study the kinetics of interfacial broadening of a polymer blend, and is validated by comparing the results to the predictions from pure Brownian dynamics and pure local dynamics calculations.


I. INTRODUCTION
Dynamical processes in soft materials are typically characterized by multiple time and length scales, an intrinsic feature that brings significant challenge for the design of efficient simulation schemes [1,2]. Although the computational power still continuously increases each year, studying macroscopic properties of polymer systems by atomistic simulations has remained infeasible. Therefore, to study processes on larger time and length scales, coarse-grained (CG) models have been developed as well as multiscale modelling schemes where different length scales are coupled concurrently. Among these simulation techniques, hybrid adaptive resolution schemes are particularly useful, since in such schemes most parts of the system are represented by a lower resolution model, to save time, while the detailed description holds only where it is required, to keep accuracy. Depending on the concrete representation schemes, the hybrid models may be categorized into hybrid particle-particle methods (HPP) [3][4][5][6][7][8][9][10], where both the fine-grained and coarsegrained description are constructed using particle-based models, and hybrid particle-continuum (or particle-field) (HPF) methods, where the fine-grained molecules are resolved by particles while the coarse-grained molecules are represented by some continuum quantities which are controlled by some continuous dynamic equations [11][12][13][14][15][16][17]. Well-known examples of HPP schemes are the hybrid QM/MM methods [18] or the hybrid atomistic-CG approaches [3][4][5][6]. HPF schemes include hybrid particledensity functional schemes [16,17] and hybrid schemes linking molecular dynamics and fluctuating hydrodynamics (MD/FH) [11,12,14,19]. However, so far, the HPF schemes have either been restricted to static properties -such as, e.g., the Monte Carlo scheme proposed by us in previous work [16] -or to fluids of small molecules -such as, e.g., the MD/FH schemes where molecular fluids are connected to Navier-Stokes fluids via buffer zones [11,14,15,25,26,28].
In the present paper, the HPF will refer to adaptive resolution models where materials are represented at different levels -particle or field -in different spatially separated domains, as explained above. However, we should note that in the literature, the term "hybrid particle-field modeling" is also used for a variety of other, intrinsically different models: Originally, the idea of hybrid particle-field modeling was exploited for the study of two-component systems with one component described by a particle-based model while the other uses a field model. For example, in the late 90s, Miao et al [20] and Saphiannikova et al [21] proposed a mixed scheme to study the conformational properties of polymer brushes under shear. They utilized the continuous Brinkman equation (a variant of the Navier-Stokes equation) to describe the motion of the solvent, and Brownian dynamics (BD) to simulate polymer brushes. Similar ideas were later adopted by Sides et al [22], with the difference that they focussed on a colloid-polymer composite, and treated polymers as continuous fields and colloids as Brownian particles. The main feature of these hybrid approaches is that different resolution schemes are chosen for specified components at the very beginning, and these representations do not change during the simulation.
Ganesan and coworkers [23,24] and Milano and Kawakatsu [29][30][31] introduced another type of singleresolution "hybrid particle-field modeling", where the fields are essentially used to define the interactions between particles. The whole simulation is controlled by a particle-based BD or MD algorithm, while the fields, e.g., densities, appear as intermediate variables, i.e., they are extracted from all the particle positions and then used for the calculation of interaction forces. The introduction of these fields can be viewed as a way to efficiently deal with nonbonded molecular interactions. A similar approach was first proposed by Laradji et al [32,33], and also characterizes the "single chain in mean field theory" [34,35].
In contrast, the MD/FD approach first proposed by de Fabritiis et al [11,12] is a "true" HPF model with adaptive resolution, which which combines particle-based MD simulations and continuum-based fluctuating hydrodynamic models dynamically. In this model the stochastic equations are integrated by numerical methods of computational fluid dynamics. Mass and momentum conservations are achieved through flux balance restriction. Later, this HPF scheme was combined with an HPP (atomistic/CG) model leading to a triple-scale ap-proach [13]. Similar MD/FD schemes using other Navier-Stokes solvers have been proposed, e.g., methods based on smoothed dissipative particle dynamics [25,26] or on multiparticle collision dynamics [28].
Recently, we have developed a different HPF method which is specifically designed for large, complex molecules, i.e., polymers. It couples a particle region to a field region where the thermodynamics is described by polymer density functional theory [16,17]. The simulation algorithm uses a combination of Monte Carlo (MC) moves and field relaxation steps. We have shown how to construct such hybrid models systematically by fieldtheoretic methods. The derivation is based on the fact that the particle representation and field representation originate from the same partition function obtained from the same Edwards Hamiltonian [36][37][38][39]. Thus in principle, the particle description, the field description, and the derived hybrid description are almost equivalent and should give the same macroscopic thermodynamic properties. This equivalence has been verified within statistical errors through the study of polymer conformation properties in a solution [17] and effective interactions between two polymer grafted colloids immersed in a diblock copolymer melt [16]. However, so far the method cannot be used to study dynamical processes.
In this paper, we extend our previous static model and move one step forward to construct a dynamic hybrid HPF scheme for the simulation of systems with diffusive dynamics. The scheme ensures local mass conservation, but it neglects inertia effects and thus does not account for momentum conservation. Therefore, it is suited for studying relaxational dynamics and diffusion behavior in systems where hydrodynamic interactions are not important. More specifically, particles move according to Brownian dynamics (BD), and fields are propagated according to the equations of local dynamic density functional theory (LD) [52][53][54]. Polymers can switch between different levels of resolution on the fly. These "resolution switches" are controlled by a chemical potential like function, i.e., a tuning function, rather than by a force or potential interpolation as in other adaptive resolution models. The switch process is constructed such that local mass conservation is guaranteed. We illustrate the method at the example of a study of interface broadening in an A/B polymer blend.
The remainder of the paper is organized as follows: In Sec II, we describe our methodology at the example of a polymer solution with implicit solvent. The three components of our hybrid scheme, the BD method (II A), the LD method (II B) and identity switches (II C) are discussed separately. In Sec III, we present simulation results. We first test the constructed hybrid scheme using the polymer solution and then utilize it to study the interdiffusion process of a polymer blend. We close with a brief summary in Sec IV.

II. HYBRID PARTICLE-FIELD SCHEME: METHODOLOGY
In this section we describe in detail the construction of the HPF model. For simplicity, we focus on a solution of homopolymer Gaussian chains. The extension to more complicated polymers or copolymers and to mixtures etc. is straightforward. We consider a system containing n t Rouse polymers of length N confined in a cubic box of volume V = L x · L y · L z with periodic boundary conditions. In the following and throughout the paper, we give all lengths in units of the radius of gyration of a free Gaussian chain where v ex is the excluded volume parameter, andφ(r) ≡ 1 ρ0 mj δ(r − R mj ) is the normalized microscopic density with respect to the total mean monomer density ρ 0 = n t N/V .
In the HPF scheme, the chains are first partitioned into "particle type chains" (p-chains) and "field type chains" (f-chains), depending on the state of virtual internal degrees of freedom as described below in Sec II C. The f-chains are then turned into fluctuating density fields by a field-theoretic transformation and a saddle-point integration [16,17,40,41]. This results in an effective free energy functional of the form (2) for the particle-field system, which depends on the positions {R mj } of p-chain monomers -the index m refers to the molecule, and the index j to the monomer in the chain -and on the local densities φ f of f-chain monomers. The terms F 0,p and F 0,f include the contributions of the intrachain interactions and the chain entropy of p-chains and f-chains, respectively, and F I describes the bonded interactions, which is derived from (1) and given by in a homopolymer solution. Hereφ p (r) of course refers to the normalized local density of p-chains only.
For Gaussian chains, F 0,p can be written as [46] The sum m runs over the p-chains and R mj is the position of the jth bead of the mth chain. Finally, the intrachain contribution of f-chains to the free energy, F f 0 , takes the general form where ω f (r) are conjugate fields defined such that hypothetical external potentials W = 1 N ω f would generate exactly the f-density configuration φ f (r) in a reference system of noninteracting f-chains, and Q f is the partition function of a single chain subject to such potentials W f . For Gaussian chains, it is given by We note that in this formalism, φ f is the only free field, and ω f must be determined from φ f . An explicit relation between φ f and ω f can be obtained through the following procedure: We define an end-integrated single chain propagator q f (r, s), which satisfies the modified diffusion equation [41,42] where s is a normalized chain contour variable ranging from 0 to 1. The propagator of q f is then calculated numerically with initial condition q f (r, 0) = 1 and periodic boundary conditions, e.g., using a pseudospectral method. Knowing the propagator q f , one can calculate the single chain partition function via Q f = 1 V dr q f (r, 1), and the density field via [41,42] φ The functional inversion to get ω f from φ f (density targeting problem) is done numerically using iteration techniques. Several iteration schemes are available; for more details, see, e.g., [43][44][45].
In the dynamic HPF scheme, the motion of the particle-resolved p-chains is described by BD dynamics, whereas the density profiles of the coarse-grained f-chains evolve according to a dynamic density functional, the LD scheme. The switches between p-chains and f-chains are controlled by a predefined tuning function. Apart from the switches, the system is propagated by numerical integration of the BD equations of motion (p-chains) and the LD equation (f-chains). More concretely the evolution of the system in each integration step is composed of three independent steps. First, all the beads' positions of p-chains are updated according to the BD equations, second, density profiles are evolved according to the LD equations, and third, a finite number of identity switches are carried out controlled by a Monte Carlo (MC) criterion. An advantage of completely decoupling the updating and the switching process is that the hybrid algorithm becomes more flexible with respect to incorporating other dynamic models, and can be extended to other polymer models more easily. In the following we describe the BD dynamics, the LD method, and the switching algorithm respectively.

A. Particle Propagation: Brownian Dynamics
The motion of particle beads (in p-chains) is driven by a conservative force derived from the Hamiltonian and a random force originating from the thermal fluctuations, i.e., . The random force f mj is Gaussian distributed with zero mean and variance f mjI (t)f nkJ (t ′ ) = δ mn δ jk δ IJ where I,J denote the Cartesian components. The derivative of Hamiltonian over particle positions can be performed directly using the chain rule and properties of the delta function. The result can be formally written as [47] whereω p can be viewed as the conjugate potential toφ p , and is obtained asω in our homopolymer solution. Therefore the nonbonded force acting on one bead is evaluated through the derivative of densities with respect to the bead's position.
In practice, all densities are discretized, and they are extracted by a "particle-to-mesh" assignment function. To do so, we uniformly divide the simulation box into n x · n y · n z cells, and define all densities on the vertexes (mesh points) of these cells. Several choices of the assignment function [48] are possible. The zeroth order scheme is the so called the nearest-grid scheme [49], which assigns a bead to its nearest mesh point. Here we adopt a linear order scheme, which assigns fractions of a bead to its eight nearest mesh points, the so-called cloud-in-cell (CIC) scheme [50,51]. In the CIC scheme, the fraction assigned to a given vortex is proportional to the volume of a rectangular whose diagonal is the line connecting the partition position and the mesh point on the opposite side of the cell. The assignment function completely determines the finite-differential form of the derivative of ω p . The nonbonded force acting on a bead is finally expressed as the linear interpolation ofω p defined on the bead's eight nearest mesh points with weights calculated according to the distance between the particle and the mesh point. As the explicit form of the assignment function and the nonbonded force are extensively discussed elsewhere, we just refer to Refs. [47,48].

B. Field Propagation: Local Dynamic Density Functional
The density fields φ f of f-chains are propagated according to a dynamic density functional. The starting point is the continuity equation for the monomer densities ρ 0 ∂φ f ∂t = −∇j, where the current j is driven by the thermodynamic force ∇µ f with µ f := δF δφ f [44]. Here we use a local dynamics scheme, where the thermodynamic force and the flux are linearly related, j = −D c N ρ 0 φ f ∇µ f , leading to the LD equation [52][53][54]: with In the present work, the LD equation is integrated numerically using the explicit Euler scheme, and the derivatives are evaluated in Fourier space using fast Fourier transform (FFT) [55]. To perform the numerical calculations, we discretize all the continuum quantities and define them only on the mesh points for convenience. Thus the space decomposition used in the LD dynamics is the same as that used in the BD simulation when calculating the densities.

C. Chain Identity Switches
Originally, all n t chains in the system are identical. For the HPF modeling, we need to first partition these chains into p-chains and f-chains. There are several ways to do the partitioning. The strategy we adopt here is based on the introduction of an additional "spin" variable τ mj , which is attached to each bead as an identifier [16,17]. The spin τ mj can take a value either 0 or 1. A chain is identified as an f-chain if all the spins on this chain are set zero, otherwise it is a p-chain. We note that this partioning method leads to the asymmetry of p-chains and f-chains in the spin space, where p-chains have a larger spin entropy. To control the value of τ mj , we further introduce a tuning function ∆µ(r), which plays the role of the spin's conjugate potential. The tuning function couples to τ mj , and thus determines the distributions of p-chain and f-chains. Therefore, p-chains and f-chains are treated as two chemically different species, and ∆µ acts similar to a chemical potential difference, which however varies locally. From a technical point of view, we thus work in a semi-grand canonical ensemble.
In the present study, the tuning function is imposed externally and must be specified beforehand. Some general considerations can help to construct an appropriate form for ∆µ. For example, it should be chosen such that p-chains appear only in specific regions (e.g., near boundaries), where a detailed description is required, while in the remaining large region, chains are represented by coarse-grained fields. Based on such considerations, an optimal form of ∆µ may be constructed that it is in compatible with the geometry of the localized region as well as the whole system.
The chain identity switches are driven by the tuning function, and to sample the spin variable, we adopt a Monte Carlo type scheme [56,57]. It is motivated by the following consideration: We require that a spin variable τ on a bead at position r assumes the value τ = 1 with a probability w 1 roughly given by [16,17] while the probability for τ = 0 is w 0 = 1 − w 1 . In the absence of other dynamical processes (BD or LD moves), this probability distribution can be generated exactly with the following algorithm: 1. We generate a random number in the interval [0, n p · M + n f , where M is an integer used to increase the probability of choosing a p-chain, which is introduced to at least partly account for the large spin entropy of p-chains. If the integer generated is smaller than n p · M , we choose a p-chain, otherwise we manipulate an f-chain.
2. Suppose that the mth p-chain is chosen (with probability M np·M+n f ), we then randomly single out a bead j in this chain (the probability that the ith bead being chosen 1/N ), and subsequently we flip the bead's identity.
If τ mj = 0, we attempt to change it to τ mj = 1. According to the definition of the p-chains, in this case, the p-chain is still a p-chain after the switch. We write the transition probability of such a switch as where W pp 01 denotes the acceptance probability, which will be specified below.
If τ mj = 1, we try to flip it to τ mj = 0. The chain type (p-chain or f-chain) after the switch depends on the spin values of all other beads. We have to distinguish between two different cases.
In the first case, at least one other bead (say the lth bead) carries a spin τ ml = 1. After the switch, the chain is then still a p-chain. The transition probability is where W pp 10 is again an acceptance probability. In the second case, we have τ ml = 0 for all l = j, and the switch will change this p-chain to an fchain. The transition probability for such a switch is written as Here, P pf j and the acceptance probability W pf j may depend on the conformation {R mj } of the chain m.
3. If an f-chain has been picked in step 1 (with probability n f np·M+n f ), we attempt to switch it to a pchain. To do so we first choose a bead j (with probability 1/N ) and a location r according to a probability H j (r), and then generate a trial bead at the position r. In case the move is accepted, this bead becomes the j-th bead of a p-chain and replaces the j-th bead of an f-chain. Hence the probability of picking that f-bead is H j (r)/ρ 0 φ f,j (r)∆V , where ∆V is the volume of a cell, and ρ 0 φ f,j the density of j-th monomers of f-chains. Then we construct a Gaussian chain with jth bead fixed at r. The a priori probability P priori to generate a given set of coordinates {R k } is where the normalization factor Z 0 = ( N 4π ) 3(N −1)/2 corresponds to te partition function of an ideal noninteracting reference chain of length N . The transition probability from an f-chain to this p-chain is thus expressed as with the acceptance probability W fp j , which again may depend on the conformation of the new pchain. 4. The trial move is accepted with probability W pp 01 , W pp 10 , W pf j , or W fp j , respectively. To ensure local mass conservation, this is done in a manner that the total densityφ p (r) + φ f (r) is preserved everywhere. Hence, if a switch from the mth p-chain to the mth f-chain is accepted, we remove its previous contribution to the density, sayφ m fromφ p , and add it to φ f . This also entails the calculation of the new ω f which gives φ f +φ m . Similarly, if the switch from an f-chain to a p-chain is accepted, we remove the density distribution of the newly generated pchain from φ f and update ω f accordingly.
If only switch moves are made, the target spin distribution (13) can be obtained by choosing the acceptance probabilities W pp 01 , W pp 10 , W pf mf , and W fp j such that they satisfy the detailed balance condition. In the case of p-p switches, this condition simply reads w 0 P pp 01 = w 1 P pp 10 , and it can be implemented with the simple Metropolis criterion W pp 01 = min(1, e −∆µ ), W pp 10 = min(1, e ∆µ ).
In the case of p-f or f-p switches, one must account for the entropy gain associated with the replacement of an actual p-chain m with density distributionφ m by a fchain density distribution φ m . The corresponding free energy difference is given by , where φ f denotes the density distribution of f-chains without the chain m. Thus the detailed balance condition reads which can be achieved with the Metropolis form ).
With this choice of acceptance probabilities, one would obtain the target spin distribution (13) in a simulation that only includes spin switches, i.e., chain identity switches. However, the expressions (21) are quite complicated and involve time consuming calculations of ∆F 0,f and φ f,j . On the other hand, an exact implementation of a target spin distribution is not truly necessary, since the spins are auxiliary variables and have no physical meaning. Moreover the BD and LD moves described in Secs II A and II B distort the spin distribution anyway, since they do not satisfy detailed balance with respect to (13). Therefore, we have chosen to simplify the expressions (21) by making the following approximations: First, we replace ∆F 0,f ≈ − ln(Z 0 ). Second, we neglect the fact that ρ 0 φ f,j depends on j, i.e., we approximate ρ 0 φ f,j ≈ 1 N ρ 0 φ j , and choose H j (r) ≡ H(r) = ρ 0 φ f (r) in step 3. With these simplifications, the acceptance probability (21) takes the simple form in analogy to Eq. (19).

A. Chain partitioning in a polymer solution
We first test the partitioning of p-chains and f-chains driven by the tuning function in a homogeneous polymer solution with periodic boundary conditions. We choose the total number of chains n t = 10 4 , and chain length N = 20 for all chains. Initially, we set n p = n f = 5000, and the configurations of p-chains are generated randomly as Gaussian coils, while the field density is set to be homogeneous. The size of the box is L x = L y = 2, and L z = 16 with discretization n x = n y = 8, and n z = 64, i.e., we have uniform cubic cells with volume 0.125. The average number of beads in each cell is about 50, which should be large enough to generate smooth densities. The excluded volume parameter is set v ex = 1.
To integrate the dynamic equations numerically, we choose a multiple time step approach with time step ∆t BD = 5 × 10 −4 t 0 /N for the BD scheme, and ∆t LD = 10 −4 t 0 for the LD equation. During each evolution step, we thus update f-chains one time and p-chains 4 times, i.e. we have ∆t ≡ ∆t LD = 4∆t BD . In each simulation step, after updating the p-chains and f-chains, we attempt to switch beads S = 500 times, and about 10 chains are successfully switched. The parameter used to enhance the probability of finding a p-chain is set to M = 20. With these parameters, the density profiles for p-chains and f-chains reach equilibrium within about half the relaxation time t 0 of a single chain.
In principle, the tuning function can be chosen at will. In our test, we assume that ∆µ depends only on z and has an almost steplike profile, are parameters determining the shape of ∆µ. Figure 1 shows the density profiles of p-chains, f-chains, and their sum after equilibration of the system. As expected, p-chains aggregate near the borders at z = −8 and z = 8, while f-chains mainly appear in the middle of the system. P-chains and f-chains are separated by interfaces with width about 2.5R g . The interface region can become even narrower if the tuning function is chosen sharper. However, due to the finite extension of chains, it must have a minimum width around R g . In principle, the total density should be independent of ∆µ, i.e., the summation of p-chain and f-chain density should be homogeneous. Indeed, in Figure (1), the total density in the p-chain region and f-chain region is almost constant apart from small statistical fluctuations. However, in the interface region, a small dip appears. Such dips are often observed in the interfacial regions of mixed resolution models [16,58,59]. In our case, it can partly be related to a discretization effect, as it becomes weaker if the number of grid points is increased. If necessary, it could be removed by introducing additional potentials in the interface region, following Ref. [58].
The efficiency of partitioning of course depends on the frequency of chain switch attempts. This is demonstrated in Figure (1) (b) and (c), where the number of attempts has been reduced from the "default" value S = 500 (Figure 1 a) to S = 50 and S = 20 (Figure 1 b) and c). The interdiffusion of p-and f-chains competes with the chain switching, and therefore, the interfacial region broadens. This also reduces the density dip.

B. Interface broadening in an A/B polymer blend
Next we test our hybrid model on a more complex problem, the interdiffusion of miscible A/B polymers at an A/B interface. For comparison, we also run pure particle-based BD simulations (considered as the reference system) and pure LD calculations for the same system.
We consider systems containing n A = 5000 A polymers and n B = 5000 B polymers with the same chain length N = 20. The melt is assumed to be compressible with κ = 10, and A/B monomers are taken to be fully compatible, i.e., the incompatibility parameter χ (the Flory Huggins parameter) is set to χ = 0. Thus the nonbonded interactions are described by the Hamiltonian where n t = n A + n B is the total number of polymers in the system. This replaces (1) in Sec II, and Eq. (3) changes accordingly.
All simulations are implemented in three dimensional space with size L x = L y = 2, and L z = 24 decomposed into n x · n y · n z = 8 · 8 · 96 uniform cells. To integrate the dynamical equations, we choose ∆t BD = 5 · 10 −4 t 0 /N , and ∆t LD = 10 −4 as in Sec III A. Similarly, we also set M = 20 and attempt 500 p-f switches are tried in each evolution step. Final statistical quantities are obtained by averaging over the results from 64 separate independent runs. Initially, we set the number of p-chains and f-chains for A and B polymers to equal values, i.e., n Ap = n Af = n Bp = n Bf = 2500. The initial configuration for p-chains are generated from additional MC simulations leading to homogeneous density profiles along the x and y directions and sharp interfaces along z direction. Therefore, in the following, we only consider the density profiles along the z direction. For simplicity, the initial density profiles of f-chains are set such that φ Af (t = 0) =φ Ap (t = 0), and φ Bf (t = 0) =φ Bp (t = 0). At t = 0, we thus have a strongly phase separated system with an A rich region near the boundaries (located at z = −12, and z = 12), and a B rich region in the middle (See Fig. (4) for the initial density profile of A polymers φ A (t = 0). The profiles of φ B (t = 0) are not shown).
We are interested in the A/B interface regions and want to keep the chain configuration properties there. Thus we need to resolve polymers as particles in the interface regions, while we can represent them by fields elsewhere. This can be done by choosing a tuning function for both A and B polymers such that it has large values in the interface regions and small values elsewhere. Specifically, we use where µ a = −10, µ b = 0, η = 10, ∆ = 1 are parameters controlling the shape of ∆µ. The tuning function ∆µ is plotted in Figure 2. successful partitioning of chains according to the tuning function. These one dimensional densities are obtained by averaging over the corresponding three dimensional densities in the x and y directions. In the early stage, t ∼ 0.1t 0 , the shape of the initial densities does not change very much, and most of the A polymers are still located near boundaries, while the B polymers are confined to a slab in the middle of the system. As expected, p-chains of both A and B polymers are mainly found in the interface regions (around L z = −6 and L z = 6) with width about 2.5R g , and f-chains are distributed elsewhere. We find that that at t ∼ 0.1t 0 , the number of p-chains is roughly n Ap ∼ n Bp ∼ 1000, while the number of f-chains is n Ap ∼ n Bp ∼ 4000, i.e., about 20% of A polymers (and B polymers) are resolved by BD particles. Small dips can also be seen in the total density profiles similar to that in polymer solutions. With the presently chosen parameters, f-chains are not fully replaced by pchains in the A/B interface regions, i.e., the maximum value of φ Ap and φ Bp does not reach the bulk value. One could easily adjust the parameters such that a wide slab in the interfacial region is filled by p-chains only. However, here we are interested in the interplay of particles and fields and its influence on the interdiffusion, which is why we keep the p-domains very narrow and f-chains can intrude. Next we investigate the interfacial broadening process at early stages. This is done by monitoring the evolution of polymer densities. Since A polymer and B polymers are symmetric, we just focus on the properties of A polymers. Figure 4 shows the evolution of the A polymer density φ A (t) =φ Ap (t)+ φ Af (t). Starting from the initial sharp density profile, φ A (t) gradually becomes diffuse. At late times (not hown), φ A (t) approaches a uniform distribution with magnitude about 0.5.
To characterize the interfacial broadening more quantitatively, we introduce a quantity W measuring the interfacial width, and monitor W as a function of time. The "interfacial width" W is simply defined as the inverse of the maximum slope of the density profile. We do not renormalize it to the "bulk densities", as the latter are changing with time and are thus not well defined. Figure 5 plots the evolution of width W (t) obtained with our hybrid method, and compares with with data from pure particle-based BD simulations and from pure fieldbased LD simulations. The initial densities in the pure LD calculations is set the same as that in the pure BD simulations. Except for early times (t < 0.1t 0 ), where deviations are observed due to the finite interfacial width effect, the curves for the hybrid model, the pure BD simulations and pure LD theory are rather well described by a scaling relation [60][61][62] W ∼ t a with exponent a ∼ 0.45.
The simulation data obtained with the hybrid model are in excellent agreement with those obtained from the fully fine-grained BD model. In contrast, the data from the LD simulations show deviations at early times. These presumably reflect the importance of density fluctuations in the interface region, which are neglected in the coarsegrained LD approach. This problem can be avoided by resorting to the hybrid description.

IV. SUMMARY AND REMARKS
In this work, we have developed a hybrid scheme that dynamically couples a particle-based description of polymers to a field-based description in an adaptive resolution sense, where polymers can switch between different resolutions on the fly. The resolution switch is controlled by a predefined tuning function through a Monte Carlo criterion. In order to capture the dynamic behavior, we use Brownian dynamics at the particle level, and a local density functional theory at the field level. The hybrid scheme was explained in detail at the example of a polymer solution, and then tested in simulations of a more complex problem, the interface broadening in a polymer A/B blend.
In the applications presented here, the tuning function was kept constant throughout the simulation. The geometric partitioning into particle and field domains was thus pre-defined prior to the simulation. However, the partitioning has no physical meaning, and the tuning function can be chosen at will, which means that it can also be adjusted to the configuration during the simulation. For example, in a large-scale simulation of A/B interfaces in immiscible polymer blends, capillary undulations of the interface must be taken into account. The fine-grained p-domain should thus be able to follow the local position of the interface. This could be done, e.g., by coupling the tuning function to the local concentration difference of A and B monomers, or to concentration gradients.
We should note apart from LD, a number of other dynamic density functional theories (DDFTs) have been proposed for dynamic studies of polymer systems [44,53]. Unfortunately, in contrast to simple liquids [63][64][65][66][67], it is not possible to derive d DDFT for polymers which is based on densities only and fully equivalent to the BD description [68,69]. Therefore, DDFTs for polymers al- The dashed line is plotted to guide eyes and it indicates a power law W (t) ∼ t a with a ≃ 0.45. The size of error bar for the data obtained from hybrid simulations is similar to that from pure BD simulations (not shown). ways rely on approximations and represent different aspects of the polymer motion. The local DDFT used here (the LD method) focusses on individual monomers and disregards the effect of chain connectivity. Other DDFTs such that Debye dynamics and external potential dynamics [53] account for the chain connectivity and focus on the collective motion of monomers in a chain. In a recent work [70], we have compared several DDFT approaches with BD simulations, and shown that the LD scheme reproduces the interfacial evolution in A/B homopolymer or the kinetics of microscopic phase separation in A:B diblock copolymer melts fairly accurately. However, in other situations, the kinetic pathways of structure formation in DDFT simulations may depend on the choice of the DDFT [71,72]. Therefore, to complete the hybrid modelling scheme, we also tested models where the f-chains are propagated according to DDFT equations that account for chain connectivity, e.g., Debye dynamics. Unfortunately, we found that such schemes tend to suffer from a flow imbalance in the f-p interfacial region. As a result, the overall density is no longer constant, an accumulation of chains is observed in the p-or f-domain, and corrections have to be applied. A detailed analysis will be published elsewhere.
Another important extension in future work will be to add inertia and account for momentum conservation. At the coarse-grained particle level, this is easily done by switching from Brownian Dynamics to Molecular Dynamics, possibly in combination with a momentum-conserving thermostat that restores friction effects [73]. At the field level, one option is to use a recently developed scheme that couples DDFT with Lattice Boltzmann simulations [71,72].
Hybrid particle-field schemes can be useful in a variety of contexts: Firstly, they can be used to implement proper boundary conditions in field-based simulations [17]. This is particularly important in situations where complex surfaces or interfaces critically influence the relevant properties of a material. Examples are liquid crystal cells where surfaces define the orientation of the director, or composite materials with polymer-coated filler particles. Secondly, they can be used for particlebased simulations of polymeric systems with open boundary conditions. The main purpose of the field domain is then to provide a realistic environment for the system of interest and act as polymer reservoir. We believe that our HPF approach and suitable refinements and extensions will provide a convenient and practical tool for such applications.