1 Introduction

Continuum mechanics was formulated more than two centuries ago and since then it has been applied successfully to study various materials. One of the main tools of continuum mechanics is the Navier-Stokes (NS) equations, which are continuum partial differential equations (PDEs) formulating conservation laws for mass and momentum of a Newtonian fluid (e.g., Refs. ^[1] and ^[2]). Over a century ago, it was realized that matter consists of discrete molecules. The complex structures and interactions among a large number of molecules manifest continuum behaviors at the macroscopic level. With advances of computer technology, both NS equations and molecular trajectories can be discretized to study hydrodynamic phenomena at distinct scales, from which numerical methods such as computational fluid dynamics (CFD) and molecular dynamics (MD) have emerged. Since the former method is fast for simulating macroscopic phenomena and the latter method elaborates on microscopic details, the two methods are not meant to compete but to complement each other. With the miniaturization of micro-fluidics technology and growing interest of fluid phenomena in the biophysics context, the mesoscopic fluid regime has been receiving more and more attention. In this spatial-temporal regime, MD is computationally expensive while a classical CFD is physically inaccurate. As a consequence, several mesoscopic methods have been developed over the last three decades to exploit a proper trade-off between molecular details and computational efficiency. Among the popular mesoscopic methods are coarse-grained molecular dynamics (CG-MD) methods^[3], lattice-Boltzmann method^[4-5], dissipative particle dynamics (DPD)^[6-7], and multiple particle colliding dynamics/stochastic rotation dynamics^[8-9]. CG-MD is a generic name, under which there are many techniques to generate a CG potential of mean force, such as reverse Monte Carlo^[10], iterative Boltzmann inversion^[11], force-matching method^[12-13], relative entropy method^[14], and Mori-Zwanzig projection^[15-16], to name but a few.

In this work, we discuss hydrodynamics of DPD exclusively, but the insights and perspective may be relevant to other mesoscopic methods as well. DPD was invented by Hoogerbrugge and Koelman^[6] 25 years ago, who postulated pairwise dissipative and random forces between a set of Lagrangian particles to comply with the isotropy and Galilean invariance of a simple fluid. Later, Espa nol and Warren^[7] added in the pairwise conservative, reformulated the updating algorithm as stochastic differential equations (SDEs), and derived the corresponding Fokker-Planck equations (FPEs). Moreover, they obtained the celebrated fluctuation-dissipation theorem (FDT) for a DPD system so that the Gibbs canonical distribution is the steady state solution of the FPEs. Thereafter, DPD has found a large variety of applications in simple and complex fluids^[17-24]. Meanwhile, researchers have been trying to establish rigorous foundations for DPD from both bottom-up coarse-graining^[25-28] and top-down discretization^[29-30]. Given an arbitrary input settings of DPD, qualitative results on transport coefficients and equation of state of hydrodynamics can also be derived from both kinetic theory^[31] and projection techniques^[32]. It is well known from both theories and simulations that DPD reproduces hydrodynamic behavior of NS equations at sufficiently large scales^[32-35]. Furthermore, DPD considers also thermal fluctuations, which emerge due to finite thermodynamic volume and are hallmarks of mesoscopic phenomena. As the mesoscopic scale is quite vaguely defined and one usually does not explicitly define the applicability range of DPD. If this arbitrariness is of concern, one may employ smoothed dissipative particle dynamics (SDPD)^[29], which discretizes the NS equations and considers thermal fluctuations consistently^[36]. Correspondingly, SDPD represents a discrete version of the Landau-Lifshiftz-Navier-Stokes (LLNS) equations^[1], and the DPD method may be considered as a "poor" version of SDPD method from this top-down perspective. In this work, however, we wish to discuss about the possibility that DPD may represent generalized hydrodynamics^{[33, 37-38]} to certain extent, which is not covered by SDPD.

In Section 2, we first introduce the classical DPD method. To have a specific reference system for discussion, we employ a toy model system of star polymer described by a Hamiltonian MD^[27-28]. To reproduce static and dynamic properties of MD, we devise three versions of DPD by two techniques: the first technique is to run DPD simulations for empirical parameters with trials and errors until we have desired properties generated by DPD; the second technique is to coarse-grain MD trajectories for DPD force fields using the Mori-Zwanzig projection^{[15-16, 25-28, 39]}. With the second technique, we consider two levels of degrees of freedom in the resultant DPD. In Section 3, we introduce the current correlation functions (CCFs), which are significant in time-dependent processes. Meanwhile, we also revisit solutions of CCFs at sufficiently large scales derived from linearized NS equations. In Section 4, we compute CCFs from MD and DPD simulations, and compare the results with hydrodynamic solutions. The aim is to identify to what extent the results of DPD agree with those of MD and continuum hydrodynamics. In the end, we summarize our findings in Section 5.

2 DPD

DPD was originally introduced by Hoogerbrugge and Koelman^[6]. Here, we present the reformulated version of Espa nol and Warren^[7], who introduced the conservative force into the original DPD,

(1)

where R_I and P_I are the position and momentum of particle I, respectively. The mass for each particle is taken as an identical constant M. For isotropy and Galilean invariance^[6], the forces are postulated to be pairwise and depend only on relative position R_IJ=R_I-R_J and relative velocity V_IJ=V_I -V_J=(P_I -P_J)/M of two particles. To preserve linear momentum, forces are antisymmetric and satisfy Newton's third law, that is, F_IJ=-F_JI. To preserve angular momentum, forces are always along the radial direction of two particles e_IJ=R_IJ/R_IJ, with R_IJ=|R_IJ| being their relative distance. Relevant letters are capitalized for a DPD system to distinguish from its underlying MD system appearing later in Subsection 2.1. Particle index I ranges from 1 to the total number of particles N.

According to the general coarse-graining theory or the Mori-Zwanzig projection formalism^{[15-16, 39-41]}, a generalized Langevin equation (GLE) results from coarse-graining, and three types of forces emerge at the mesoscopic level of description. Therefore, the total force on each particle in DPD is decomposed into three parts which are postulated to be^[6-7]

(2)

where α, γ, and σ are strengths of individual forces. Weighting functions w_C, w_D, and w_R are isotropic and depend only on R_IJ. Beyond a cut-off radius R_c, values of w_C, w_D, and w_R vanish. ξ_IJ=ξ_JI is a Gaussian white noise, that is,

(3)

where δ_IJ is the Kronecker delta function and δ(t-t') is the Dirac delta function. If we substitute the forces of Eq. (2) into Eq. (1), we can write the GLE in a mathematically well-defined form of SDE as

(4)

where dW_IJ=dW_JI are independent increments of the Wiener process, and Itȏ calculus is assumed. Therefore,

(5)

where dW_IJ is an infinitesimal of order 1/2 in time^{[7, 18]}. Therefore, the total force in Eq. (1) should be rewritten as .

The FPEs corresponding to the SDEs of DPD can be formulated to acquire a steady state solution of the canonical form^[7]. It turns out that if the following relations

(6)

are satisfied, a DPD system has an invariant canonical distribution. Equation (6) is the celebrated FDT of DPD, which was first derived by Español and Warren^[7].

The central task for a successful DPD simulation is to design a proper setting of R_c, α, γ, w_C(R), and w_D(R), which may be probed by various techniques. To anchor our discussion on a quantitative analysis, we consider a polymer melt system modeled by MD as reference. This idealized MD system is then taken as the "first principle" in this work. Therefore, any version of DPD should refer to the results of MD to evaluate its performance. At first, we shall present the details of the MD system. This is followed by elaborations of three versions of DPD resulting from two types of coarse-graining techniques.

2.1 MD

We consider an MD system of a polymer melt at temperature k_BT=1.0, with k_B being Boltzmann's constant. The Hamiltonian defining the phase space trajectories is the sum of the kinetic energy and potential energy of the system

(7)

where Γ^p and Γ^r are momentum and position spaces, respectively. N_t is the total number of microscopic/atomic particles. Each polymer is a molecule with chains of beads connected by short springs. The polymer is isotropic and has a star shape with a center bead connecting N_a arms with N_b monomers per arm. Therefore, each polymer contains N _c=N_a× N_b+1 microscopic particles. For simplicity, we only consider N_a=10, N_b=1, and N_c=11 in this work.

The potential energy is the sum of pairwise potentials, each of which has two components as V(r)=V_WCA(r)+V_FENE(r). The non-bonded interaction is modeled by a truncated Lennard-Jones (LJ) potential, also named as Weeks-Chandler-Andersen (WCA) potential^[42],

(8)

for which r_c=2^1/6σ so that only the repulsive part of the LJ potential remains. The bonded interaction is modeled by a finitely extensible nonlinear elastic (FENE) potential^[43],

(9)

with bond strength and bond reference length r₀=1.5σ. We further set the length scale σ=1, the energy scale , the mass of each bead m=1, and therefore, the time scale . We take the simulation box to be cubic with length L≈ 30.184 1, and fill N=1 000 molecules into the box. Therefore, the total number of microscopic particles is N_t=1 000× N_c=11 000, leading to a number density ρ=N_t/L³≈0.4.

The trajectories of microscopic particles in phase-space are determined by the Hamiltonian/Newtonian mechanics as

(10)

In simulations, the trajectories are integrated by the measure-preserving Verlet method^[44] with time step δt=10^-3τ. For a canonical ensemble of MD simulations, we apply the virial equation^[37] to calculate the pressure of the system at equilibrium, which measures as P=0.191. The diffusion coefficient for the center of mass (CoM) for each polymer measures as D=0.119, by either the time derivative of its mean-squared displacement (MSD) or the integral of its velocity autocorrelation function (VACF). The kinematic viscosity of the MD system measures as ν=0.965, with either a small-perturbed reversible Poiseuille flow^{[28, 45]} or transversal current correlation function (TCCF)^[34-35]. These values have also been reported in a previous work^[28]. Results from different measuring techniques agree with each other.

2.2 Empirical DPD

To achieve certain properties of a target fluid by the classical DPD method, input parameters R_c, α, γ (or σ), and weight functions w_C(R) and w_D(R) (or w_R(R)) are calibrated empirically, that is, one performs simulations by trials and error until satisfactory properties are generated. Here, we adopt typical weight functions as suggested^{[18, 20]}

(11)

DPD simulations are performed in a cubic box with the same size as that of MD. Then, the number density of DPD particles is reduced to be ρ/N_c≈0.036 4. Here, we fix the cut-off radius as R_c=3.4 as suggested in Subsection 2.3 before calibrating other input parameters. After a series of simulations by trials and error^[28], we set α=13.5 to reproduce the pressure of MD. We further determine γ=32.0 and s=0.5 to reproduce the diffusion coefficient and viscosity of MD. This set of parameters leads to a DPD system with P=0.194±0.04, D=0.119± 0.002 and ν=1.012 3±0.06¹, closely resembling the properties delivered by the MD simulations. All quantities in DPD are reported also in LJ units of MD.

¹There is a typo in Ref. ^[28] where it was reported as ν=0.444 on Table 3.

Note that the empirical input parameters of DPD may also be probed in a more systematical approach with a stochastic optimization process^[46-47]. Since radial distribution function (RDF) and VACF for particles in the empirical DPD were not set as the optimized target in the beginning for deducing the input parameters, it is not surprising that they deviate from those of MD^[28]. About the effects of arbitrariness from R_c and other input parameters on RDF and VACF of the DPD system, we refer to another thorough study^[46]. These later subtleties, however, are not of particular concern of the current study and we shall not discuss them further.

2.3 Mori-Zwanzig guided DPD

Using an elegant projection technique named as Mori-Zwanzig formalism in non-equilibrium statistical mechanics^{[15-16, 39]}, one may derive closed-form dynamic equations for a set of CG variables, which are functions of microscopic phase space Γ^{[25, 27]}. Here, we give a brief introduction following a derivation given in Ref. ^[25].

2.3.1 Mori-Zwanzig formalism

In the case of a polymer melt, the CG variables of interest are quantities associated with the CoM of each polymer molecule such as its position and translational momentum,

(12)

(13)

We further define the phase-space density for the CG variables as

(14)

where the CG phase-space coordinate is defined as Γ_s={R₁, R₂, ..., R_N, P₁, P₂, ..., P_N} and are the corresponding field variables. The probability distribution function at equilibrium for the microscopic phase space Γ is taken as the canonical ensemble, that is, with β=1/(k_BT) and Z^-1 being the normalizing partition function. Therefore, the probability distribution function at equilibrium for the CG phase space Γ_s is

(15)

We define configuration and momentum distribution functions as

(16)

(17)

where and are the position and momentum spaces of CG field variables, respectively, and the integrals are always for N_t dimensions. Then,

To construct a Hilbert space for dynamical variables, we define the scalar product as the correlation function in the canonical ensemble at equilibrium,

(18)

where dynamical variables are considered to be real so that B^*(Γ) =B(Γ). If we utilize the set of phase-space density functions as the basis in the Hilbert space, and further require the completeness condition^{[25, 40]}

(19)

then any dynamic variable A(Γ) can be decomposed into two parts by two complementary projection operators as

(20)

Since we expect to be expanded by the basis in the Hilbert space, the projection operator is naturally defined as

(21)

Therefore, the definition for is , or

Applying the definition of projection operators, we may decompose the evolution equation , and thereafter obtain . The last term by definition is the evolution equation for the CG momentum and is expressed as^[25]

(22)

with the random forces defined as

(23)

(24)

(25)

Here, δF_I is the fluctuating force, which equals the difference between the instantaneous total force F_I and mean force on particle I. To verify that the first term on the right-hand side of Eq. (22) is the mean force, we have

(26)

which is the mean force. We also note the essential difference between δF_I(t) and , where the former evolves with time by the usual microscopic propagator , whereas the latter evolves with time by the orthogonal propagator and stays always orthogonal to the projected subspace. By integrating Eq. (22) in time, we have information on P_I, and therefore can determine the position trajectory dR_I=P_I dt for each CG particle as well.

Although Eq. (22) is insightful and exact, it is still generally difficult to apply for a practical simulation, as it is an integro-differential equation involving many-body effects, with non-Markovian terms and projected dynamics. Firstly, we assume that the mean force is pairwise decomposable and many-body terms are neglected,

(27)

Therefore, the pairwise term can be calculated as a function of distance between CoMs in an ensemble of MD simulations. Secondly, for the convolution (memory) term, we make a Markovian approximation, i.e.,

(28)

(29)

This approximation is justified, as the CoMs are slow variables compared with the MD^[28].

The last assumption is to approximate the projected dynamics by the real dynamics, that is, . To calculate Γμν, if we simply took the integral for the correlation of real dynamics up-to ∞ as

(30)

we would run into the plateau problem^[48], i.e., the integral vanishes with ^{[16, 49]}. However, if we take the integral up-to τ as

(31)

where an inequality holds. Thus, we do have a plateau region and achieve a definite good approximation of Γ_IJ^{[16, 49]}. In practice, we apply Eq. (31) with τ dynamically determined so that the integral always reaches the maximum value of all possible τ. In this way, we calculate Γ_IJ as a function of distance between CoMs in the ensemble of MD simulations.

If the rotation of the CG particle is of interest, the angular momentum can be defined for each polymer molecule as

(32)

where the averaged moment of inertial is I_I=6.55 for the star polymer and Ω_I is the angular velocity. Accordingly, the evolution equation of L_I is also available.

2.3.2 Parameters of DPD

Following the Mori-Zwanzig formalism described above, we obtain DPD parameters as α=795.69, γ=γ_||=146.18, and weight functions as w_C(R)=(1+4R/R_c)(1-R/R_c)⁴ and w_D(R)=w_D^||(R)=(1+3R/R_c)(1-R/R_c)³ with R_c=3.32. The form and coefficients of the polynomials for the weight functions are obtained by fitting the numerical data generated by MD simulations^[28]. For reference, we term this model as Mori-Zwanzig dissipative particle dynamics or MZ-DPD for brevity. Given γ and w_D(R), we may obtain σ and w_R(R) from the FDT of Eq. (6). Finally, we generate Gaussian white noises for the pairwise random force of MZ-DPD to perform an ensemble of DPD simulations.

It is remarkable that in the MD simulations of polymer melt, the instantaneous forces between different polymer molecules are not along the radial direction e_IJ of their CoMs, but in all directions. Since the conservative force depends only on relative positions, its average vanishes in all directions but the direction along e_IJ due to symmetry. However, the dissipative force depends not only on relative positions, but also on relative velocities. Therefore, we need to introduce the transversal component in addition to the radial component between two molecules to represent the extra degree of freedom in the MD system^[50]. Meanwhile, rotational degrees of freedom for each DPD particle also need to be considered so that both linear and angular momenta of the system are conserved at the CG level. The extra transversal and rotational dissipative forces have already been hypothesized in a phenomenological way in a previous work^[50]. For the full dynamics of the CG variables, we term the model as Mori-Zwanzig full DPD, or MZ-FDPD for brevity, and the equation of motion reads^[28]

(33)

(34)

where T_I is the torque for particle I. Here, γ_||, w_D^||, w_R^||, and σ_|| are for the dissipative and random forces, respectively, along the radial direction e_IJ, while and are for the dissipative and random forces, respectively, on the plane perpendicular to e_IJ. dW_IJ is a matrix of independent Wiener process with dW_IJ^A being its antisymmetric part. Therefore, based on the parameters and weight functions for MZ-DPD, we further have =110.76 and =(1+3.95R/R_c)(1-R/R_c)^3.95 with R_c=3.32. It is worth noting that MZ-DPD is just a simplified model by keeping only the forces along the radial direction e_IJ between two particles in MZ-FDPD. MZ-DPD does not consider the equation of motion of angular momentum either.

As the conservative force determines the static properties, both MZ-DPD and MZ-FDPD produce the same pressure of P=0.193, close to that of MD. However, due to lack of the extra degree of freedom in the dissipative force, MZ-DPD underestimates the kinematic viscosity as ν=0.851 and overestimates diffusion coefficient as D=0.138. In contrast, the MZ-FDPD with all necessary degrees of freedom represented on the CG level. It resembles the MD closely with D=0.12 and ν=0.954. Before further discussing about the performance of different versions of DPD, we introduce the concept of correlation functions in the next section, which is crucial in time-dependent processes.

3 CCFs

Temporal correlation functions play a similar role in time-dependent processes as the partition function in equilibrium thermodynamics. For a comprehensive discussion on time-dependent phenomena, we refer to Refs. ^[37], ^[38], and ^[49]. In particular, here, we are concerned with autocorrelation functions, i.e., correlations between the same dynamical variable at different times. More precisely, if A(t)=A[r₁(t), r₂(t), ..., r_N(t), p₁(t), p₂(t), ..., p_N(t)] is a dynamical variable of coordinates and momenta of all particles of a system, its autocorrelation function is defined as

(35)

where t=t'-t"≥ 0 and A^* is the complex conjugate of A.Angular brackets denote equilibrium ensemble average or time average. We assume that the system is ergodic so that both averages are equivalent. At equilibrium, the system is stationary so that the correlation functions are invariant under time translation. Therefore, C_AA(t', t")=C_AA(t=t'-t"). Autocorrelation functions are real functions of time^{[37, 49]}.

For a hydrodynamic variable described by the hydrodynamics equations such as mass density or momentum density, which is a continuous fields variable, we have the corresponding microscopic expression of the dynamic variable in an MD or DPD system. The microscopic definition is

(36)

where a_i is the mass or momentum of particle i. According to Eq. (36), the microscopic mass density and momentum density are defined as

(37)

(38)

where mass for each particle is a constant m and u_i is the velocity for particle i. The corresponding continuous variables are related to the microscopic variables as

(39)

(40)

To distinguish, we denote the continuous dynamic variable with an over-bar as A(r, t). The definitions of the continuous variables require a "coarse-graining" explicitly in space over a volume v. The size of v must be macroscopically small to allow for a continuum definition. Meanwhile, it must be microscopically large so that the relative fluctuations of the number of particles within v are negligible^[37]. When v is sufficiently large, Eqs. (39) and (40) imply a coarse-graining in time as well. However, an exact value of v is problem-dependent and is not definitely specified in general.

To compare the autocorrelation functions among the results of MD, DPD of different versions, and hydrodynamic equations, we should avoid the ambiguity of the same dynamic variable in physical space, which may be interpreted differently by the three distinct descriptions. Therefore, we consider the Fourier-transformed counterpart of the dynamical variable in k-space defined as

(41)

where the integral is taken over the whole volume V of the system and k is the wave vector. Since it is translation invariant at equilibrium, there is no correlation between A_k and A_k' when k≠k'^[51]. Accordingly, the corresponding autocorrelation function in k-space with the same k is defined as

(42)

In particular, when the dynamic variable is A_k(t)=j_k(t), the autocorrelation function is a 2nd rank tensor and it is named as CCFs. Due to isotropy, it is diagonal and may be further decomposed into transversal and longitudinal components as

(43)

where α, β=x, y, or z and are Cartesian components of the unit vector . Without loss of generality, we take k=(k, 0, 0)^T. Therefore, the longitudinal and transversal current correlation functions (LCCFs and TCCFs), are expressed as

(44)

More specifically, if we consider the dynamics at continuum, that is, lim k → 0, we may ignore the difference between j_k(t) and j_k(t) and consider C_T(k, t) ≈ C_T(k, t) and C_L(k, t) ≈ C_L(k, t)^[37-38]. Therefore, we shall neglect the over-bars on the dynamical variables without causing ambiguity at lim k → 0. Moreover, the continuum CCFs can also be expressed by the solutions of linearized hydrodynamic equations. We shall revisit such solutions as follows.

3.1 Hydrodynamic solutions

A dynamic variable is conserved if it satisfies a continuity equation in a differential form as

(45)

where j^A is the current associated with the variable A. Therefore, the conservation laws for mass and momentum are expressed as

(46)

(47)

where j is the mass current or momentum density introduced in Eqs. (38) and (40). Π is the momentum current or stress tensor, and for a Newtonian fluid, it is defined as

(48)

where P(r, t) is the local pressure, and η and ζ are dynamic shear and bulk viscosities, respetively. If we choose a frame of reference where the mean velocity vanishes, that is, , Eqs. (46) and (47) may be linearized by considering the variations of the hydrodynamic variables as

(49)

(50)

(51)

A variation of second order is dropped leading to one "≈" sign in j(r, t). Inserting Eq. (51) into Eqs. (46) and (47), we obtain

(52)

(53)

Furthermore, ∇P(r, t)=∇(P +δP(r, t))=∇ δP(r, t)=c_T²∇ δρ(r, t), with c_T being the isothermal sound speed. The above equations can be conveniently solved by performing a space-Fourier transform into k-space introduced in Eq. (41),

(54)

(55)

where the shear kinematic and longitudinal kinematic viscosities are introduced as

(56)

(57)

It is apparent from Eq. (54) that the evolution of density ρ _k(t) is coupled with the longitudinal component of the momentum density j_k^x(t)=k·j_k(t), but decoupled from the transversal components. Therefore, we do not consider separately, as it provides merely a redundant information as C_L(k, t). This observation is valid not only at lim k → 0, but at an arbitrary k, as Eq. (54) utilizes the mere fact of mass conservation without any other assumptions. Therefore, we may rewrite Eqs. (54) and (55) as four equations of scalar variables

(58)

(59)

(60)

(61)

Equations (60) and (61) are first-oder linear homogeneous ordinary different equations (ODEs), which are readily solved as^[52]

(62)

(63)

Equations (58) and (59) are two coupled ODEs and may be written in a matrix form as

(64)

where a and H are defined as

(65)

(66)

The general solutions to Eq. (64) are determined by the eigenvalues of H, which can be obtained by further solving det (H-λI})=0^[52]. It is easy to show that the eigenvalues are

(67)

with

(68)

where the sound attenuation coefficient is introduced as Γ_T. Here, we consider an under-damped solution, where s_T is real, that is, k < 2c_T/ν_L. In particular, we consider a continuum limit where . Therefore s_T≈ c_T and we arrive at the solutions as

(69)

(70)

The normalized temporal autocorrelation functions of the fluctuating variables in k-space are readily expressed as

(71)

(72)

where is assumed to derive C'_L(k, t). In the subsequent section, we always calculate the normalized quantities C'_T(k, t) and C'_L(k, t) and therefore omit the sign "'" without ambiguity. It is worth noting that the TCCFs and LCCFs are self-similar at continuum scales and they scale with powers of k. Therefore, self-similar results of MD or DPD at different k would indicate a continuum behavior, whereas different properties at different k resort to a non-continuum behavior.

4 Results and discussion

Here, we calculate the CCFs for the polymer melt model introduced in Section 2 from simulations by MD, empirical DPD, Mori-Zwanzig guided DPD of two versions, and solutions of the linearized hydrodynamic equations. Thereafter, we are able to compare the results from DPD simulations of different versions, against that of MD and of hydrodynamic solutions. Since we are focusing on the CCFs in Fourier-transformed k-space, we avoid the difficulty of one-to-one correspondence between hydrodynamic variables expressed by particles in MD and hydrodynamic variables expressed in the other two coarser descriptions. The input parameters for MD and DPD simulations are as given in Section 2. The wave vector is taken as k_n=(k_n, 0, 0) with k_n=2πn/L, where n is a positive integer number and L is the length of a cubic box. We will consider n=1, 2, 3, and 4, where for the first two wave numbers, a continuum behavior is apparent, while for the latter two wave numbers, a particle behavior is prevalent.

4.1 Continuum limit

We compare TCCFs of MD and continuum theory of Eq. (71), where the latter adopts ν=0.965 as in MD simulations. From Fig. 1(a), it is clear that at length scales L and L/2 (wave number k₁ and k₂ accordingly), MD simulations behave like continuum hydrodynamics and its TCCFs follow closely exponential decays, where the decay rate is νk². In continuum hydrodynamics, the determining parameter of TCCFs is the kinematic viscosity ν=0.965 (η=0.386, ρ=0.4). Therefore, the classical DPD with empirical input parameters for targeting the MD viscosity produces almost identical TCCFs as those of MD and continuum theory, as shown in Fig. 1(a). However, with the same number of degrees of freedom as that of the empirical DPD, MZ-DPD underestimates the decay rate of TCCFs, as shown in Fig. 1(b). This is not surprising given the fact that MZ-DPD underestimate the viscosity of MD, as reported already in Subsection 2.3. In contrast, with all necessary degrees of freedom on the CG level, MZ-FDPD reproduces well the viscosity of MD, and therefore also the TCCFs of MD, as indicated in Fig. 1(c).

To recapitulate, TCCFs of DPD resemble those of MD and hydrodynamic solutions at continuum scales (k₁ and k₂) with the same exponentially decaying law characterized by the kinematic viscosity. This is true as long as the version of DPD reproduces the same viscosity of MD or continuum theory, as the viscosity is the only determining factor for TCCFs in the continuum scales and TCCFs are self-similar with exponential-decay rates νk².

For the LCCFs, the situation is more complicated than that for TCCFs, as the former present themselves with both decaying and oscillating behaviors characterized by Γ_T and c_T, as shown in Eq. (72). We may fit well the hydrodynamic solution (72) to MD results, as indicated in Fig. 2(a). After the fitting, parameters Γ_T=1.17 and c_T=1.3 are determined as well. As further by-products, we determine ζ=0.42, ν_L=2.34. Since the empirical DPD is developed to target the pressure and viscosity of the MD system, it is not expected to reproduce ζ and c_T of MD well. This deficiency is reflected by the mismatch on LCCFs between DPD and MD in Fig. 2(a). In contrast, the Mori-Zwanzig guided DPD and full DPD reproduce the LCCFs of MD very well, as indicated in Figs. 2(b) and 2(c). It is also interesting to note that there is no apparent difference between the performance of MZ-DPD and MZ-FDPD on LCCFs, although MZ-FDPD has non-radial interactions between particles. This indicates that the longitudinal modes (Γ_T and c_T) are purely determined by the radial-direction interactions between particles.

We note that for the approximation made in Eq. (68) for s_T≈ c_T, we have assumed that k₁≈ 0.208 or k₂≈ 0.416 2c_T/ν_L ≈ 1.11 is valid. From the results shown above, the assumption of "" does not lead to an apparent deficiency of the hydrodynamic solutions on LCCFs.

4.2 Below continuum limit

If we further consider smaller scales k₃=6π/L and k₄=8π/L, we clearly observe that the decay of TCCF for MD is no longer strictly exponential, as indicated in Fig. 3. Here, the solutions of hydrodynamic equations are still exponential with kinematic viscosity taken as ν=0.965. Since we take MD as the first principle, its results are reliable and should be taken as reference. Therefore, the hydrodynamics solutions are incorrect at these small scales. From the development of generalized hydrodynamics^[37-38], it is not difficult to realize that we have a crude assumption for our hydrodynamic equations; the fluid is assumed to be Newtonian at all scales in Eq. (48), that is, for α ≠ β. This is a Markovian relation between the shear stress and instantaneous shear rate. For fluids (not only this star polymer system, but also simple fluids such as argon) at sufficiently small scales, there is in general a non-Markovian relation, which determines the shear stress from the history of shear rates. This is to say that the fluid is viscoelastic at small scales^{[37-38, 53]}. Therefore, if one wants to compare hydrodynamic solutions with that of MD, one needs to model this non-Markovian relation phenomenologically (with memory kernels such as exponential, and Gaussian functions) or determine the non-Markovian relation from the projection technique introduced earlier^[54-56]. Once the memory kernel is available, it is possible to obtain accurate hydrodynamic solutions given that the non-Markovian hydrodynamic equations are still solvable. This is a common practice during the development of generalized hydrodynamics and we shall not proceed further in this work.

It is remarkable that the TCCFs of the empirical DPD resemble those of MD quite well at k₃ and k₄, especially for intermediate and long time (t 2), as indicated in Fig. 3. This is surprising as the interacting laws between DPD particles are Markovian and the non-exponential behavior of TCCFs is non-Markovian in origin. It appears that certain coherent movements of DPD particles at length scales of 2π/k₃=L/3 and 2π/k₄=L/4 manifest themselves as non-Markovian effects in the TCCF and show viscoelastic effects. We need to mention that the viscoelastic effects from MD at these small scales are also due to the coherent movements of MD particles. So far the analogy, if there is any relation at all, between the TCCFs of MD and DPD at k₃ and k₄ is still not clear.

Moreover, the TCCFs of MZ-FDPD resemble those of MD even better, except for t→ 0, as indicated in Fig. 4. We note that the interacting laws between particles in MZ-FDPD are also Markovian and their differences from empirical DPD are the actual values of coefficients and weight functions. These results suggest that even with empirical Markovian DPD, we do have the possibility to adjust input settings wisely so that we recover non-Markovian/viscoelastic/generalized hydrodynamic behaviors of a reference MD system.

For the LCCF, the results can neither be recovered by the hydrodynamic solutions nor by the empirical DPD, as indicated in Fig. 5. For the deficiency from the hydrodynamic solutions, the first questionable assumption is the same as in the case of TCCF that a Markovian relation between stress and strain rate is invalid, and furthermore η and ζ should be also scale-dependent. Even if the Markovian relation was correct, the assumption about s_T≈ c_T is invalid, as it implies that k₃≈ 0.624, k₄≈0.833 2c_T/ν_L ≈ 1.11. For the disagreement from the empirical DPD simulations, it is not surprising, as we never intend to recover the LCCF when deriving the empirical interactions between the particles. The disagreements were already observed for k₁ and k₂. However, it is remarkable to note that for these small scales at k₃ and k₄, the LCCFs of MZ-FDPD represent reasonably well of that of MD, as indicated in Fig. 6, except at very small times (equivalently higher frequencies). This suggests again that it is possible to adjust input settings of DPD so that generalized hydrodynamic behavior of a reference MD system may be recovered.

5 Summary

We present three versions of DPD; the first version is empirical so that the pressure, viscosity, and diffusivity of a reference system are simultaneously well reproduced. The force fields of the second and third versions are obtained from a systematical coarse-graining procedure within the Mori-Zwanzig projection formalism by keeping the properties of the CoM in MD descriptions. The difference between 2nd and 3rd versions are that the former considers only radial-direction interactions between neighboring particles, just as the empirical/classical DPD formulations, while the latter considers in addition transversal and rotational interactions between particles. They are named MZ-DPD (Mori-Zwanzig DPD) and MZ-FDPD (Mori-Zwanzig full DPD), respectively.

To have some quantitative comparison, we consider a specific system of star-polymer melt at number density ρ=0.4 modeled by MD as a reference system, which we consider as the first principles in this work. Since the coarse-graining level is N_c=11, that is, eleven atoms are coarse-grained into one DPD particle, time scales between CG variable and MD dynamics are well separated. Therefore, we may replace the memory kernel with a Dirac delta function, and furthermore the friction forces between DPD particles are Markovian and random forces are white-colored. We also approximate the random force, which evolves by an orthogonal propagator, by the usual propagator describing the real dynamics. This is justified for the Markovian description, as long as we calculate the friction coefficients by taking the maximum value from the integral of time correlations of random force. As a result, the typical plateau problem is avoided.

Previously, we have evaluated the performance of DPD of three versions such as pressure, viscosity, diffusion, and velocity autocorrelation function of the CoM of the polymer, in comparison with results of MD in physical space^[28]. If we want to compare the performances of DPD further with hydrodynamics solutions, we would have to introduce the definition of hydrodynamic field variables for MD and DPD. This involves an arbitrary length scale as shown in Eq. (40), which is generally difficult to define. Therefore, in physical space, it is difficult to compare fairly the dynamic processes at different scales described by MD, DPD, and hydrodynamic equations. However, in the Fourier-transformed k-space, we can focus on a particular wave number or length scale at each time and compare the three descriptions in detail. For such purposes, we adopt and revisit the concept of the CCFs. From the calculations of transversal and longitudinal CCF, or TCCFs and LCCFs for brevity, we observe that at continuum scales, results of DPD of three versions resemble overall well those of MD and hydrodynamic solutions for a Newtonian fluid. One notable disagreement is that LCCFs of empirical DPD deviate from those of MD and hydrodynamic solution. This is not surprising, as the target for optimizing the DPD interactions does not include sound speed c_T or longitudinal viscosity Γ_T. Another notable disagreement is that TCCFs of MZ-DPD deviate slightly from those of MD. This is due to the underestimation of MD viscosity from MZ-DPD, which does not consider transversal or rotational interactions between neighboring particles. The empirical DPD also only considers radial-direction interactions between neighboring particles, but does not have a problem with TCCFs. This is because the friction coefficient γ and weight function w_D(R) in empirical DPD are free to tune so that it compensates the missing effects of transversal/rotational interactions on the viscosity, while the friction coefficient and weight function in MZ-DPD are completely determined by the projection procedure.

Not only for the polymer fluids of small molecules considered here, but also for simple fluids such as argon, the dynamical behavior is viscoelastic/non-Markovian at small scales^[38]. Therefore, the hydrodynamic equations for a Newtonian fluid do not describe the dynamics accurately at small scales. This is reflected by the TCCFs and LCCFs calculated at k₃ and k₄, where hydrodynamics solutions deviate significantly from those of MD. However, the MZ-DPD and MZ-FDPD are two Markovian approximations of the GLEs. It is remarkable that the Markovian interacting laws reproduce quite well the non-Markovian behavior of the TCCF and LCCF from MD. Even more remarkable is that the empirical DPD which has Markovian interaction laws as well, can also reproduce that of MD at small scales, represented by the comparison of TCCF between them. These phenomena are no coincidence and we think that the coherent movements² of DPD particles should be explored more so that even an empirical Markovian DPD may recover the generalized hydrodynamics to some extent. We report our first observations here and further investigations are on the way.

²Coherent movements represent the coordinated motion of particles at scales larger than individual particles and correlations between velocities of particles which are significant

Acknowledgements acknowledges funding support of the U. S. Army Research Laboratory with Cooperative Agreement No. W911NF-12-2-0023.