Hydrodynamic relaxations in dissipative particle dynamics

TL;DR

This study investigates relaxation phenomena in dissipative particle dynamics using fluctuating hydrodynamics, revealing temperature-dependent accuracy of hydrodynamic predictions and the dominance of sound wave propagation at high temperatures.

Contribution

It provides a detailed analysis of hydrodynamic relaxations in DPD, highlighting temperature effects and the independence of transverse dynamics from shear force contributions.

Findings

Hydrodynamic theory predicts transverse dynamics well at low temperatures.

Agreement with hydrodynamics depends on viscosity definition.

High temperatures favor sound wave propagation in longitudinal spectra.

Abstract

This paper studies the dynamics of relaxation phenomena in the standard dissipative particle dynamics (DPD) model [Groot and Warren, JCP, 107:4423 (1997)]. Using fluctuating hydrodynamics as the framework of the investigation, we focus on the collective transverse and longitudinal dynamics. It is shown that classical hydrodynamic theory predicts the transverse dynamics at relative low temperatures very well when compared to simulation data, however, the theory predictions are, on the same length scale, less accurate for higher temperatures. The agreement with hydrodynamics depends on the definition of the viscosity, and here we find that the transverse dynamics are independent of the dissipative and random shear force contributions to the stress. For high temperatures, the spectrum for the longitudinal dynamics is dominated by the Brillouin peak for large length scales and the relaxation is therefore governed by sound wave propagation and is athermal. This contrasts the results at lower temperatures and small length scale, where the thermal process is clearly present in the spectra. The Landau-Placzek ratio is lower than the classical model Lennard-Jones liquid, especially at higher temperatures. The DPD model, at least qualitatively, re-captures the underlying hydrodynamical mechanisms, and quantitative agreement is excellent at intermediate temperatures for the transverse dynamics.

Tables

Table 1. Table 1: Table of the viscosities, η 0 1 superscript subscript 𝜂 0 1 \eta_{0}^{1} and η 0 2 superscript subscript 𝜂 0 2 \eta_{0}^{2} , the self-diffusivity, D s subscript 𝐷 𝑠 D_{s} , the Schmidt number, Sc, the Maxwell relaxation time, τ M subscript 𝜏 𝑀 \tau_{M} , and instantaneous shear modulus, G ∞ subscript 𝐺 G_{\infty} . The two values for the Schmidt number are for η 0 1 / ( ρ ​ D s ) superscript subscript 𝜂 0 1 𝜌 subscript 𝐷 𝑠 \eta_{0}^{1}/(\rho D_{s}) and η 0 2 / ( ρ ​ D s ) superscript subscript 𝜂 0 2 𝜌 subscript 𝐷 𝑠 \eta_{0}^{2}/(\rho D_{s}) . The uncertainties associated with the viscosities are the standard deviation of the mean calculated from five independent simulations. There are no statistical uncertainty on the digits for D s subscript 𝐷 𝑠 D_{s} , Sc. G ∞ subscript 𝐺 G_{\infty} is calculated from the sample averages of η 0 1 superscript subscript 𝜂 0 1 \eta_{0}^{1} and τ M subscript 𝜏 𝑀 \tau_{M} with one significant decimal.

$T$	1.00		0.80		0.60		0.40		0.20		0.10
${\eta }_0^1$	0.715 $\pm$ 0.006		0.661 $\pm$ 0.008		0.673 $\pm$ 0.004		0.778 $\pm$ 0.004		1.425 $\pm$ 0.008		4.13 $\pm$ 0.03
${\eta }_0^2$	0.859 $\pm$ 0.009		0.82 $\pm$ 0.01		0.848 $\pm$ 0.005		1.00 $\pm$ 0.01		1.80 $\pm$ 0.02		4.83 $\pm$ 0.07
$D_s$	0.300		0.230		0.159		0.089		0.028		0.006
Sc	1/1		1/1		1/2		3/4		17/21		229/268
${\tau }_M$	0.075 $\pm$ 0.003		0.075 $\pm$ 0.002		0.082 $\pm$ 0.003		0.105 $\pm$ 0.004		0.20 $\pm$ 0.01		0.33 $\pm$ 0.01
$G_{\infty }$	9.5		8.8		8.2		7.4		7.1		12.5

Equations

$$\frac{\partial \rho \varphi }{\partial t}={\sigma }^{\varphi }-\nabla {\mathbf{J}}^{\varphi }-\nabla \cdot (\rho \varphi \mathbf{u}),$$

$$\frac{\partial \rho }{\partial t}$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}$$

$$\frac{\partial \rho e}{\partial t}$$

$$\rho \mathbf{u}=({\rho }_{\text{av}}+\delta \rho )\delta \mathbf{u}\approx {\rho }_{\text{av}}\delta \mathbf{u}\text{and}\rho e\approx {\rho }_{\text{av}}\delta e.$$

$${\mathbf{J}}^m$$

$$\mathbf{P}$$

$${\mathbf{J}}^e$$

$$\frac{\partial }{\partial t}\delta \rho$$

$${\rho }_{\text{av}}\frac{\partial }{\partial t}\delta \mathbf{u}$$

$${\rho }_{\text{av}}\frac{\partial }{\partial t}\delta e$$

$$\delta p={\left(\frac{\partial p}{\partial \rho }\right)}_T\delta \rho +{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\delta T=\frac{1}{{\rho }_{\text{av}}{\chi }_T}\delta \rho +\frac{{\beta }_V}{c_V}\delta e,$$

$${\sigma }_e=\frac{T{\beta }_{\text{V}}}{{\rho }_{\text{av}}}\frac{\partial \delta \rho }{\partial t}=\frac{T{\beta }_{\text{V}}}{{\rho }_{\text{av}}}\left(D{\nabla }^2\delta \rho -{\rho }_{\text{av}}\nabla \cdot \delta \mathbf{u}-\nabla \cdot \delta {\mathbf{J}}^m\right).$$

$$\tilde{f}(\mathbf{k},t)={\iiint }_{-\infty }^{\infty }f(\mathbf{r},t)e^{-i\mathbf{k}\cdot \mathbf{r}}\mathrm{d}\mathbf{r}$$

$$\frac{\partial }{\partial t}\widetilde{\delta \rho }$$

$${\rho }_{\text{av}}\frac{\partial }{\partial t}\widetilde{\delta \mathbf{u}}$$

$${\rho }_{\text{av}}\frac{\partial }{\partial t}\widetilde{\delta e}$$

$$\frac{\partial }{\partial t}{\widetilde{\delta u}}_x$$

$$\frac{\partial }{\partial t}{\widetilde{\delta u}}_z$$

$$\frac{\partial }{\partial t}\widetilde{\delta \rho }$$

$$\frac{\partial }{\partial t}{\widetilde{\delta u}}_y$$

$$\frac{\partial }{\partial t}\widetilde{\delta e}$$

$$C_{AB}(\mathbf{k},t)=\frac{1}{V}⟨A(\mathbf{k},t)B(-\mathbf{k},0)⟩,$$

$$\frac{\partial C_{uu}^{\perp }}{\partial t}=-{\nu }_0{k}^2{C}_{uu}^{\perp }$$

$$C_{uu}^{\perp }(k,t)=\frac{k_{B}T}{{\rho }_{\text{av}}}{e}^{-{\nu }_0{k}^{2}t},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{ccc}{C}_{\rho \rho }& C_{\rho u}& C_{\rho e}\\ C_{u\rho }& C_{uu}& C_{ue}\\ C_{e\rho }& C_{eu}& C_{ee}\end{array}\right]=-\left[\begin{array}{ccc}D{k}^2& i{\rho }_{\text{av}}k& 0\\ \frac{ik}{{\rho }_{\text{av}}^2{\chi }_T}& {\nu }_l{k}^2& \frac{ik{\beta }_V}{c_V{\rho }_{\text{av}}}\\ \frac{T{\beta }_{V}D{k}^2}{{\rho }_{\text{av}}^2}& \frac{iT{\beta }_{V}k}{{\rho }_{\text{av}}}& \kappa k^2\end{array}\right]\left[\begin{array}{ccc}{C}_{\rho \rho }& C_{\rho u}& C_{\rho e}\\ C_{u\rho }& C_{uu}& C_{ue}\\ C_{e\rho }& C_{eu}& C_{ee}\end{array}\right].$$

$$\{C_{\rho \rho },C_{u\rho },C_{e\rho }\},\{C_{uu},C_{\rho u},C_{eu}\},\text{and}\{C_{ee},C_{\rho e},C_{ue}\}$$

$$C_{AB}(\mathbf{k},t)=K_1{e}^{-D_T{k}^{2}t}+e^{-\Gamma k^{2}t}\left[K_2\cos (c_{s}kt)+i{K}_3\sin (c_{s}kt)\right]$$

$$D_T=\frac{\kappa }{{\chi }_T{\rho }_{\text{av}}{c}_s^2}\mathrm{and}\Gamma =\frac{1}{2}\left[\frac{\kappa }{{\chi }_T{\rho }_{\text{av}}{c}_s^2}+(D+{\nu }_l+\kappa )\right]$$

$$c_s^2=\frac{{\beta }_V^2{\chi }_{T}T-{\rho }_{\text{av}}{c}_V}{{\chi }_T{c}_V{\rho }_{\text{av}}^2}$$

$$C_{AB}^N(\mathbf{k},t)=K_{AB}{e}^{-D_T{k}^{2}t}+(1-K_{AB})e^{-\Gamma k^{2}t}\cos (c_{s}kt).$$

$$\frac{\mathrm{d}{\mathbf{r}}_i}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\mathbf{p}}_i}{\mathrm{d}t}$$

$${\mathbf{F}}_{ij}^C=a_{ij}(1-r_{ij}){\hat{\mathbf{r}}}_{ij},$$

$${\mathbf{F}}_{ij}^R=\frac{\sigma w(r_{ij}){\zeta }_{ij}}{\sqrt{\Delta t}}{\hat{\mathbf{r}}}_{ij}\text{and}{\mathbf{F}}_{ij}^D=-\frac{(\sigma w(r_{ij}){)}^2}{2mT}\left[{\hat{\mathbf{r}}}_{ij}\cdot ({\mathbf{v}}_i-{\mathbf{v}}_j)\right]{\hat{\mathbf{r}}}_{ij},$$

$$\tilde{\rho }(\mathbf{k},t)$$

Full text

I Introduction

The dissipative particle dynamics (DPD) method hoogerbrugge:epl:1992 ; espanol:epl:1995 is widely used to perform mesoscale computer simulations of, e.g., polymer solutions kong:jcp:1997 , spinodal decomposition novik:pre:2000 , fluid flows in micro- and nanopores boromand:cpc:2015 ; ranjith:mn:2015 , and cell membrane damage groot:biophys:2001 , just to name a few examples. A standard DPD simulation involves a set of point particles interacting by three different forces: a conservative, a dissipative, and a random force in such a manner that momentum is conserved. The DPD particle can be thought of as a collection of molecules moving in a coherent fashion espanol:simu:2002 . The forces are often tweaked to mimic specific fluidic systems, e.g., the particles can be connected with spring forces to simulate polymer solutions and melts; see also the review by Moeendarbary *et al. * moeendarbary:jam:2009 . Importantly, the interparticle conservative force is weak and usually without a strong repulsive core, in fact, the conservative force is not necessary in order to obtain hydrodynamic behavior marsh:pre:1997 ; ripoll:jcp:2001 .

In the DPD model by Groot and Warren groot:jcp:1997 , the conservative force is linear with respect to the distance between the two point masses. This model is simple and very appealing; however, it yields an unrealistic equation of state which is quadratic in density groot:jcp:1997 . Also, the dissipative force depends only on the position and velocity differences of the two interacting particles and neglects shear forces espanol:jcp:2017 . Nevertheless, the parameter space for this model is quite large and the physical interpretation of the parameters is not always straightforward. For example, the particle density can be chosen as a free parameter for a given system, and from this choice the conservative force parameter can be estimated using the compressibility groot:jcp:1997 . Interestingly, this so-called adaptive parameter approach leads to a decreasing viscosity for decreasing temperature boromand:cpc:2015 , which characterizes a gas mcquarrie:book:1976 . This gaseous behavior is also manifested by a Schmidt number of order unity groot:jcp:1997 , where the Schmidt number is defined as the ratio between the kinematic viscosity and the diffusion coefficient. Bocquet and Charlaix bocquet:csr:2010 conjectured that classical hydrodynamics is valid for wavevectors $k$ fulfilling $k<\sqrt{2\pi\rho/\eta_{0}\tau_{s}}$, where $\rho$ is the density, $\eta_{0}$ the shear viscosity and $\tau_{s}$ is the relaxation time given by the shear stress relaxation hansen:langmuir:2015 . From this criterion one can see that in the low density limit (low Schmidt number) the classical hydrodynamic theory will break down even at large length scales as the viscosity and relaxation time are only functions of temperaure in this limit.

The hydrodynamic properties for the DPD technique have been thouroughly investigated in the past, see for example Refs. marsh:pre:1997, ; ripoll:jcp:2001, . However, as the DPD model is widely used by the simulation community keaveny:jcp:2005 ; backer:jcp:2005 ; reddy:pf:2009 ; boromand:cpc:2015 at low Schmidt number, we believe it is important to investigate the properties of the model by Groot and Warren under conditions where the Schmidt number varies from unity to higher values typically characterizing liquids like the model Lennard-Jones liquid used in classical molecular dynamics.

It has been noted by several authors espanol:pre:1995 ; marsh:pre:1997 , that the energy is not conserved in the standard DPD model, and that it cannot be applied to study systems characterized by a sustained temperature gradient on the macroscopic time scales. However, the model does feature fast energy relaxations and, as also concluded by Marsh et al. marsh:pre:1997 , it can indeed be applied to investigate these relaxations. We wish to include this here as it will provide valuable insight into the underlying mechanisms of the DPD method in general.

We base our investigation on Onsager’s regression hypothesis, which states that the regression of microscopically induced fluctuations in equilibrium follows the macroscopic laws of small non-equilibrium disturbances onsager:physrev:1931 , i.e., thermally induced perturbations relax according to hydrodynamics. Typically, these (fast) relaxations do not refer to hydrodynamic quantities like density and momentum directly, but instead to the decay of the associated correlation functions kadanoff:annphys:1963 , as predicted by hydrodynamic theory. We derive these correlation functions from basic fluctuating hydrodynamics theory as this may not be known to the reader; also, we present it in a slightly different form (albeit equivalent) from that of standard texts alley:pra:1983 ; hansen:book:2006 ; boon:book:1991 . To make the study manageable, we focus on a limited part of the parameter space of the standard DPD model.

II The hydrodynamic relaxation functions

In general, one can write the balance equation for any hydrodynamic quantity per unit mass $\phi=\phi(\mathbf{r},t)$ at position $\mathbf{r}$ and time $t$ as degroot:book:1984

[TABLE]

where $\mathbf{u}$ is the streaming velocity, $\sigma^{\phi}$ the production term, and $\mathbf{J}^{\phi}$ the flux of $\phi$. In the case $\sigma^{\phi}=0$ the quantity is locally conserved. The hydrodynamic quantities we study here are the mass density, $\rho=\rho(\mathbf{r},t)$, the streaming velocity, $\mathbf{u}=\mathbf{u}(\mathbf{r},t)$, and the excess kinetic energy per unit mass, $e=e(\mathbf{r},t)$; the latter quantity is defined as the difference between the local and average kinetic energy per unit mass, $me(\mathbf{r},t)=E_{\text{kin}}(\mathbf{r},t)-\frac{3}{2}k_{B}T$, where $m$ is the particle mass. Based on the microscopic hydrodynamic operator formalism hansen:langmuir:2015 one can derive the following the balance equations on the form of Eq. (1) in the absence of any external driving forces

[TABLE]

where $\mathbf{J}^{m}$ is the mass flux tensor due to density gradients, $\mathbf{P}$ is the pressure tensor, and $\mathbf{J}^{e}$ the excess kinetic energy flux tensor. Importantly, the excess kinetic energy per unit mass, $e(\mathbf{r},t)$, is not a conserved quantity; hence, a production term $\sigma^{e}$ appears in Eq. (2c). Furthermore, for the mass balance equation, Eq. (2a), we have decomposed the mass flux into two parts; one due to thermal motion, $\mathbf{J}^{m}$, and one due to the fluid advective motion, $\rho\mathbf{u}$.

The three quantities can be written as the sum of the constant average part and the fluctuating part, i.e., $\rho=\rho_{\text{av}}+\delta\rho$, $\mathbf{u}=\delta\mathbf{u}=(\delta u_{x},\delta u_{y},\delta u_{z})$, and $e=\delta e$ since the averages of the streaming velocity and excess kinetic energy are zero. To first order in the fluctuations we have

[TABLE]

Using the framework of fluctuating hydrodynamics zarate:book:2006 , we now introduce the linear constitutive relations with stochastic forcing

[TABLE]

where $D$ is the mass flux diffusivity coefficient, $p$ is the normal pressure, $\eta_{v}$ and $\eta_{0}$ the bulk and shear viscosities, $\lambda$ the heat conductivity, $c_{V}$ the specific heat per unit mass at constant volume, and $\stackrel{{\scriptstyle os}}{{(\boldsymbol{\nabla}\mathbf{u})}}$ is the trace-less symmetric part of the strain rate tensor.

Equations (4b) and (4c) are just the constitutive relation for a Newtonian fluid and Fourier’s law of conduction mcquarrie:book:1976 with added stochastic forcing. However, as we cannot in general ignore cross-correlation effects on small time and length scale, it is noted that $D$ is not the self-diffusion coefficient maginn:jpc:1993 . Since the mass density and excess kinetic energy are scalars, that is of the same parity, both fluxes in Eqs. (4a) and (4c) can depend on the gradients of $\rho$ and $e$ according to Courier’s principle degroot:book:1984 . Here we follow Alley and Alder alley:pra:1983 and model the cross coupling through the production term $\sigma_{e}$ and the pressure $p_{eq}$.

In equilibrium the stochastic forcing term has a zero average zarate:book:2006 and is uncorrelated with the hydrodynamic quantities, e.g., $\langle\delta\mathbf{J}^{m}(\mathbf{r},t)\delta\mathbf{u}(\mathbf{r^{\prime}},t^{\prime})\rangle=\mathbf{0}$. Substituting Eqs. (3) and (4) into Eq. (2), we arrive at the stochastic dynamics. To first order in the fluctuations this is

[TABLE]

since the advective terms are of second order. More advanced stochastic descriptions have been developed in order to, for example, include elastic properties of the fluid tabak:jam:2015 ; wang:jsc:2016 . For local thermodynamic equilibrium, the pressure fluctuations can be written as hansen:book:2006

[TABLE]

where $\chi_{T}=-1/V\left(\partial V/\partial p\right)_{T}$ is the isothermal compressibility, $\beta_{\text{V}}=\left(\partial p/\partial T\right)_{\rho}$ is the thermal pressure coefficient, and $\delta e=c_{V}\delta T$. The production term for the excess kinetic energy is given by Alley and Alder alley:pra:1983

[TABLE]

Defining the Fourier transform as

[TABLE]

and then substituting Eqs. (6) and (7) into Eq. (5) gives, in Fourier space,

[TABLE]

If one makes a particularly simple choice for the wavevector, then the dynamics can be decomposed into transverse (normal) and longitudinal (parallel) dynamics with respect to this wavevector. For example, if we select $\mathbf{k}=(0,k,0)$, then from Eq. (9) the transverse dynamics is given by the streaming velocity components $\widetilde{\delta u}_{x}$ and $\widetilde{\delta u}_{z}$ via

[TABLE]

where $\nu_{0}=\eta_{0}/\rho_{\text{av}}$ is the kinematic viscosity. We will use both the dynamic viscosity, $\eta_{0}$, and kinematic viscosity, $\nu_{0}$, whenever one is more convenient than the other. As expected, Eqs. (10a) and (10b) are identical with respect to the dynamics and that the transverse dynamics are independent of the energy and density fluctuations. The longitudinal dynamics are given by

[TABLE]

where $\nu_{l}=(\eta_{v}+4\eta_{0}/3)/\rho_{\text{av}}$ is the longitudinal kinematic viscosity and $\kappa=\lambda/(c_{V}\rho_{\text{av}})$.

As mentioned above, one usually does not study the fluctuating quantities directly, but rather the associated correlation functions. To this end we define the equilibrium time-correlation function between quantities $A$ and $B$ as

[TABLE]

where $V$ is the system volume. Thus, multiplying Eqs. (10a) with $\widetilde{\delta u}_{x}(-k,0)$ and taking the ensemble average over initial conditions leads to

[TABLE]

for the transverse relaxation. Here $C_{uu}^{\perp}=\langle\widetilde{\delta u}_{x}(k,t)\widetilde{\delta u}_{x}(-k,0)\rangle/V$ is the transverse velocity autocorrelation function, and we have used that the stochastic forcing term is uncorrelated with the fluctuating quantities. The solution to Eq. (13) is

[TABLE]

where the initial value $C_{uu}^{\perp}(k,0)=k_{B}T/\rho_{\text{av}}$ is found from equipartition mcquarrie:book:1976 .

From Eq. (11) one can form nine coupled correlation functions for the longitudinal dynamics. For example, dynamic equations for $C_{\rho\rho},C_{\rho u},C_{\rho e}$ are formed by multiplying Eq. (11a) with $\delta\widetilde{\rho}(-\mathbf{k},0),\delta\widetilde{u}(-\mathbf{k},0)$, and $\delta\widetilde{e}(-\mathbf{k},0)$, respectively, and taking the ensemble average. In matrix notation, using the definition in Eq. (12) yields the following coupled linear differential equation system

[TABLE]

The coefficient matrix is referred to as the hydrodynamic matrix hansen:book:2006 . By performing the matrix multiplication in Eq. (15) it is seen that the longitudinal dynamics can be divided into three sets of co-dependent correlation functions, for example, $\dot{C}_{\rho\rho}=A_{1}(C_{\rho\rho},C_{u\rho})$, $\dot{C}_{u\rho}=A_{2}(C_{\rho\rho},C_{u\rho},C_{e\rho})$, and $\dot{C}_{e\rho}=A_{3}(C_{\rho\rho},C_{u\rho},C_{e\rho})$, where $A_{1},A_{2}$ and $A_{3}$ are linear functions represented by the hydrodynamic matrix. The three sets are written as triplets

[TABLE]

and each set of coupled differential equations can be solved from the hydrodynamic matrix. Up to second order in the wavevector, the solution for any of the nine correlation functions has the form

[TABLE]

where

[TABLE]

are the thermal diffusivity and sound attenuation, respectively, and $c_{s}$ defined as

[TABLE]

is the adiabatic speed of sound. The three integrating factors $K_{1},K_{2},$ and $K_{3}$ are found from the initial conditions and are, in fact, not independent. Now, $C_{AB}$ is either a real or purely imaginary valued function, which means that if $K_{3}=0$ then in general $K_{2}\neq 0$ and $K_{1}\neq 0$  while if $K_{3}\neq 0$ then $K_{2}=K_{1}=0$. In the case where $C_{AB}$ is real, the normalized correlation function is written in the form

[TABLE]

Thus, the longitudinal dynamics are governed by three fundamental processes with frequencies $D_{T}k^{2},\Gamma k^{2}$, and $c_{s}k$. From Eq. (18), one sees that $D_{T}$ pertains to the thermal processes and that the sound attenuation $\Gamma$ dampens the wave propagation with speed $c_{s}$; the magnitude of this damping is governed by all three diffusive processes, i.e., by $D,\nu_{l}$, and $\kappa$. Equations (14) and (20) form the framework for this hydrodynamics study.

III Simulation methodology

The standard DPD model by Groot and Warren is composed of a single type of point particle. The particle position, $\mathbf{r}_{i}$, and momentum, $\mathbf{p}_{i}$, follow Newton’s equation of motion,

[TABLE]

The total force, $\mathbf{F}_{i}$, is composed of the conservative force, $\mathbf{F}_{i}^{C}$, due to the interaction between the particles, a random force, $\mathbf{F}_{i}^{R},$ simulating the coarse graining of many degrees of freedom, and a dissipative force, $\mathbf{F}_{i}^{D}$, removing the viscous heating generated from the random force. Thus $\mathbf{F}_{i}=\mathbf{F}_{i}^{C}+\mathbf{F}_{i}^{R}+\mathbf{F}_{i}^{D}$. As it is common practise, we use reduced units such that the characteristic mass and length scales are set to unity. Also, temperature, $T$, is in units of $k_{B}/\epsilon$, where $\epsilon$ is the characteristic energy scale. In reduced units the conservative force is

[TABLE]

where $a_{ij}$ is a parameter that quantifies the repulsion between particles $i$ and $j$, $\mathbf{r}_{ij}$ is the vector of separation $\mathbf{r}_{i}-\mathbf{r}_{j}$, $r_{ij}=|\mathbf{r}_{ij}|$, $\hat{\mathbf{r}}_{ij}=\mathbf{r}_{ij}/r_{ij}$. Here we use $a_{ij}=25$ and the interactions are ignored when $r_{ij}>1=r_{c}$. Following Groot and Warren groot:jcp:1997 , the random and dissipative forces are

[TABLE]

where $\sigma$ is the random force amplitude, $\zeta_{ij}$ is a uniformly distributed random number with zero mean and unit variance, $w(r_{ij})$ is a weighing function given by $w(r_{ij})=1-r_{ij}$, $\mathbf{v}_{i}$ the velocity of particle $i$, and $\Delta t=0.02$ is the time step used in the integrator. In all simulations the amplitude $\sigma$ is set to 3.0. The equations of motion are integrated using the standard velocity Verlet algorithm by Groot and Warren groot:jcp:1997 . The system size is 5000 particles at density $\rho_{\text{av}}=3.0$, and temperatures (in reduced units) in the range $0.1\leq T\leq 1.0$ are simulated.

Español and Serrano espanol:pre:1999 studied the DPD model in terms of dimensionless parameters, namely, friction, $\Omega=\sigma^{2}r_{c}/(6v_{T}k_{B}Tm)$ where $v_{T}=\sqrt{k_{B}T/m}$, an overlap parameter, $s=r_{c}\rho^{1/3}$, and system length scale, $\mu=L_{\mathrm{box}}/r_{c}$. For relatively large friction and overlap the particle dynamics are affected by the surrounding fluid, that is, one would expect strong collective hydrodynamics. On the other hand, for low friction and small overlap the dynamics are characterized by single particle properties described by what Español and Serrano call kinetic theory espanol:pre:1999 . In the simulations carried out here, we only vary the temperature giving $1.5\leq\Omega\leq 14.7$, $s\approx 1.4$ and $\mu\approx 6.9$, and we span both the kinetic (high $T$) and hydrodynamic regime (low $T$).

During the simulations, all ten correlation functions are evaluated from the microscopic definition of the hydrodynamic variables, which to first order in fluctuations are

[TABLE]

The viscosity at zero wavevector and frequency is also evaluated. Recently, based on generic projection methods espanol:ptrs:2002 ; ernst:epl:2006 Jung and Schmid jung:jcp:2016 argued that the correct Green-Kubo integral is

[TABLE]

where the double index $\alpha\beta$ runs over the $xy$, $xz$, and $yz$ components of the pressure tensor; superscript 2 on $\eta$ distinguishes it from a viscosity defined by Groot and Warren groot:jcp:1997 and used below. $P_{\alpha\beta}^{C}$ are the three off-diagonal elements of the Irving-Kirkwood pressure tensor irving:jcp:1950

[TABLE]

and $P_{\alpha\beta}^{D}$ and $P_{\alpha\beta}^{R}$ are the dissipative and random off-diagonal components of the tensors

[TABLE]

Other authors have evaluated the viscosity based on the Irving-Kirkwood pressure only

[TABLE]

We will compare the predictions from the hydrodynamic theory using both definitions, Eqs. (25) and (28). The complex viscosity is calculated from the Irving-Kirkwood pressure tensor, i.e.,

[TABLE]

Finally, the self-diffusivity coefficient, $D_{s}$, is evaluated from the Green-Kubo integral of the single particle velocity autocorrelation function. We find that this leads to lower statistical uncertainties compared to evaluating $D_{s}$ using the particle mean-square displacements.

In a few cases, the dynamics of the DPD model is compared to a liquid-phase Lennard-Jones system at the state-point $(\rho,T)=(0.85,1.121)$ in units of $\sigma^{3}$ and $k_{B}/\epsilon$. The Lennard-Jones particles interact through the standard shifted 12-6 potential lennard-jones:prsla:1924 using a cut-off distance at $r/\sigma=2.5$. The system size is $N=1000$, and the equations of motion are integrated using the leap-frog method frenkel:book:1996 . To control the temperature, the Nose-Hoover thermostat nose:molphys:1984 ; hoover:pra:1985 is applied. The dynamic properties are calculated as explained above.

IV Results and Discussion

It is informative to study the fluid structure for the different state points investigated. Figure 1(a) plots the radial distribution functions for three state points, namely, $T=1.00$, 0.40, and 0.10; recall the density is always $\rho_{\text{av}}=3.0$. The structure can be compared to the corresponding transport properties in Table 1. First, one sees that the Schmidt number Sc$=\nu_{0}/D_{s}$ is around 1 for $T>0.6$ and that the viscosity decreases for decreasing temperature in the range $0.8\leq T\leq 1.0$, which is the well-known gas-like behavior mcquarrie:book:1976 and also reported by Boromand et al. boromand:cpc:2015 . In agreement with this, the radial distribution function shows very little fluid structure in this temperature region.

At the lowest temperature $T=0.10$, there is a clear fluid structure and the Schmidt number is of order 102. There are no indications that the system is crystaline for this temperature; for example, the mean square displacement does not feature any long time plateau, indicating no caging of the particles, and a fluidic diffusive behavior is observed after a short time, see Fig. 1 (b). For reference, the Lennard-Jones liquid state point is characterized by Sc $\approx$ 50. It is interesting that for $T=0.40$ a clear fluid structure is also absent in agreement with a Schmidt number of unity and a viscosity of $\eta_{0}^{1}=0.70\pm 0.01$ and $\eta_{0}^{2}=0.90\pm 0.01$ close to that of $T=1.0$.

To study the mechanical properties further we evaluate the shear modulus $G^{*}=G^{\prime}+^{\prime}iG^{\prime\prime}=i\omega\eta^{*}$; the loss modulus is plotted in Fig. 2 for $T=1.00,0.20$ and 0.10. Data are compared to a single-element Maxwell model

[TABLE]

where the Maxwell relaxation time, $\tau_{M}$, is found from the peak frequency in the data and using amplitude $G_{0}$ as fitting parameter. The instantaneous shear modulus (infinite-frequency complex shear modulus), $G_{\infty}$, can then be found from the relation $\eta_{0}/\tau_{M}=G_{\infty}$. Both $\tau_{M}$ and $G_{\infty}$ are listed in Table 1.

From Fig. 2 (a) it is seen that for $T=1.00$ and $\omega<20$ the DPD model is Maxwellian, or equivalently, that the shear relaxation follows a simple exponential decay for $t>\pi/10$. For $T=0.10$ the single-element Maxwell model breaks down at around $\omega=0.4$. As the temperature decreases $\tau_{M}$ increases, thus, the shear relaxation slows down as expected. We also test for time-temperature superposition (TTS) in Fig. 2(b). Here the frequency is scaled by a factor $a_{T}=\eta_{0}(T)/\eta_{0}(T_{\text{ref}})$ and the magnitude of $G^{*}$ by $b_{T}=T_{\text{ref}}/T$ ferry:book:1980 , where the reference temperature is $T_{\text{ref}}=1.00$. TTS applies for sufficiently low frequencies, but fails around $\omega\approx 1/\tau_{M}$. The shift factor $a_{T}$ increases by a factor of $\sim 6$ as temperature decreases by an order of magnitude.

Next, we turn to the non-zero wavevector regime. We plot in Fig. 3 (a) and (b) the transverse velocity autocorrelation for different wavevectors at temperatures $T=1.00$ and 0.40.

It is clearly seen that the hydrodynamical theory, Eq. (14), predicts the transverse relaxation dynamics very well in the low vector regime using the Irving-Kirkwood definition of the pressure tensor. Applying the Jung-Schmid definition gives too fast a relaxation, which indicates that this particular dynamical mode is not dependent on the random and dissipative shear forces. More quantitatively: the theory predicts the half-life as $t_{1/2}=\ln(2)/(\nu_{0}k^{2})$, i.e. for $T=0.4$ we have $t_{1/2}=9.5$ using $\eta_{0}^{1}=0.78$ and $t_{1/2}=7.4$ using $\eta_{0}^{2}=1.00$. This can be compared to the simulation result $t_{1/2}=9.5$. For very short times the theory fails to predict the relaxation; this is to be expected as the viscosity is in general both frequency and wavevector dependent, hence, for sufficiently short times the time dependence of the viscosity is important.

Interestingly, the agreement is less satisfactory for $T$=1.0; here the Irving-Kirkwood definition yields $t_{1/2}$ = 10.4 versus the simulation result $t_{1/2}=11.2$. In Fig. 4 we plot the mean square deviation

[TABLE]

where $C_{uu}^{\perp}(k,t_{i})$ is the predictions from the theory, and $C_{uu,k,t_{i}}^{\perp}$ simulation data. To avoid this parameter being affected by the large noise-to-signal ratio at very long times, we only sum over the $N_{s}$ times with data points $C_{uu}^{\perp}(k,t)/C_{uu}^{\perp}(k,0)\geq 0.1$.

Clearly, the mininum deviation is found within the temperature region $0.3\leq T\leq 0.7$. For higher temperatures the agreement is not as satisfactory; here we approach the kinetic regime as defined by Español and Serrano espanol:pre:1999 , that is, low friction and overlap parameters mentioned above. For low temperatures one observes a quite large deviation, especially pronounced for larger wavevectors. This, we argue, is due to the large characteristic frequency, $\omega=\nu_{0}k^{2}$, which is outside the classical hydrodynamic regime. For $T=0.1$ this hydrodynamic regime is never reached because of the limitations on the wavevector $k_{\mathrm{min}}=2\pi/L_{\mathrm{box}}$.

Fourier-Laplace transformation of Eq. (14) leads to

[TABLE]

which gives a peak in the imaginary part of the spectrum at $\omega_{\text{peak}}=\nu_{0}k^{2}$. This peak frequency found from the simulations is plotted in Fig. 5 (a) for $T=1.00$ and $T=0.40$ together with the hydrodynamic predictions. For low wavevectors, the peak frequency follows the predictions: $\omega_{\text{peak}}$ is proportional to $k^{2}$ and the relaxation is governed by the diffusion of momentum. The prediction fails for larger wavevectors; at lower temperature the deviation is significant for relatively lower wavevectors compared to high temperature. Again, we attribute this to the large characteristic frequency at low temperature and large wavevector.

The frequency and wavevector dependent shear viscosity can be defined by re-arranging Eq. (32),

[TABLE]

This result can also be derived from first principles by including the position and time dependence of the transport coefficient in Eqs. (4). In the zero frequency limit we have the viscosity kernel $\widetilde{\eta}(k)=k_{B}T/k^{2}\widehat{C}_{uu}^{\perp}(\mathbf{k},0)$. Figure 5 (b) shows this viscosity kernel at zero frequency for $T=1.00,0.40$ and $T=0.20$. The zero wavevector viscosity is also indicated using $\eta_{0}^{1}$ from Table 1. It is interesting to see that for $k$ less than unity, the wavevector-dependent viscosity reaches $\eta_{0}$, i.e., the local Newtonian law of viscosity holds for $k<1.0$. This is observed (in appropriate reduced units) for many different fluids hansen:langmuir:2015 . We also note that Ripoll et al. ripoll:jcp:2001 studied the kernel for zero conservative force.

Rather than approaching the deviation between theory and simulation through wavevector dependent transport coefficients, one can generalize the stochastic forcing and assume ${\delta\mathbf{J}}^{m}$, ${\delta\mathbf{P}}$, and ${\delta\mathbf{J}}^{e}$ to be correlated with hydrodynamic quantities. In this case the transverse dynamics are governed by the equation

[TABLE]

where

[TABLE]

Applying a Fourier-Laplace transform gives the correlation between forcing and the transverse velocity in terms of wavevector and frequency as

[TABLE]

Because the theoretical predictions are relatively large for higher temperatures the contribution from $\varepsilon$ is larger for all wavevectors compared to the intermediate temperatures.

We now turn to the longitudinal relaxation dynamics and focus first on the density autocorrelation function, $C_{\rho\rho}$. It is worth noting that this is related to the coherent intermediate scattering function, $F(\mathbf{k},t)$, by $C_{\rho\rho}(\mathbf{k},t)=\rho_{\text{av}}F(\mathbf{k},t)$. The density autocorrelation function is a real-valued function and, hence, it relaxes according to Eq. (20); it is plotted in Fig. 6 (a) for $T=1.00$ at wavevectors $k=0.53,2.12$ and 5.30.

The dashed line is the best fit of Eq. (20) to data using $K_{\rho\rho},D_{T},\Gamma$ and $c_{s}$ as fitting parameters. The damped oscillations predicted from hydrodynamics are evident, indicating sound waves that are dampened by the sound attenuation coefficient, $\Gamma$. The existence of the thermal process is less clear. To study this in more detail we investigate the corresponding spectra. The Fourier-Laplace transform of Eq. (20) is

[TABLE]

Again, note the two different contributions to the relaxation. The real part of $\widehat{C}_{\rho\rho}^{N}$ is symmetric about $\omega=0$ and we therefore only discuss the behaviour for $\omega\geq 0$. The first term relates to the thermal process and gives rise to the Rayleigh peak at $\omega=0$; this process is only present at low frequencies and the half-width of the Rayleigh peak is $D_{T}k^{2}$. The second term has a peak at frequency $\omega_{\mathrm{peak}}=c_{s}k$; the maximum is identified as the Brillouin peak and has half width $2\Gamma k^{2}$. Inspired by Hansen and McDonald hansen:book:2006 , the peaks and their widths are illustrated in Fig. 6 (c). Figures 6 (b) and (c) show the real part of the spectrum of the density autocorrelation function for wavevectors $k=0.53$ and $k=5.30$, the highest wavevector studied. Using the fitted values obtained in Fig. 6 (a), we plot the predicted spectra together with the transformed data. The agreement is not perfect as the local minimum predicted by the theory (at $\omega\approx 4$ for $k=5.3$) is not found in the spectrum of the data. Fitting to Eq. (37) did not improve this. For $T=1.0$ and relatively small wavevectors, Fig. 6 (b), the thermal process is almost completely absent and the relaxation is athermal. However, for large wavevector, Fig. 6 (c), the process is indeed observed in the spectrum.

The ratio of the two processes is quantified from the Landau-Placzek ratio cummins:jcp:1966 , that is, the ratio between the Rayleigh and Brillouin integral regions, or intensities, $I_{R}/2I_{B}=\gamma-1$, where $\gamma$ itself is the ratio between the heat capacities at constant pressure and volume, $\gamma=C_{P}/C_{v}$. The integral regions are also illustrated in Fig. 6 (c). In Fig. 7 (a) we plot the dispersion relation for $\gamma$ for different temperatures.

It is clear that the thermal process intensity increases as we decrease temperature and wavelength. For reference, the Lennard-Jones liquid features $1.6\leq\gamma\leq 2.6$ for $0.46\leq k\leq 5.9$, see also Bryk et al. bryk:jcp:2010 . In this region, the Lennard-Jones system also shows a clear de Gennes narrowing hansen:book:2006 ; we have not observed this narrowing for the wavevectors and temperatures studied here. From Fig. 6 (c) it is seen that the frequencies of two processes overlap indicating that the processes are coupled; this coupling is only present on relatively small length scales, that is, for typical simulation setups these two processes are decoupled and, furthermore, the thermal process only accounts for a small fraction of the hydrodynamic relaxation. However, for $T=0.10$, the coupling is relatively large even on longer length scales and may affect the response considerably.

The dispersion relation for the Brilluion peak frequency, $\omega_{\text{peak}}$, is plotted in Fig. 7 (b); it is seen that the oscillatory frequencies are roughly the same for the different temperatures at sufficiently small wavevector, which means that the speed of sound is to a good approximation independent of temperature on these length scales. For larger wavevectors the discrepancy between $T=0.10$ and $T>0.10$ is pronounced; the underlying mechanical reason for this is not well known, but likely due to the different local liquid structure on these small scales, see for example Ref. trachenko:rpp:2016, , but also the coupling of the longitudinal processes can be important. It is worth noting that the maxima seen in Fig. 7 (b) is also observed in the Lennard-Jones liquid. From the simulation data we cannot conclude if the DPD model features positive or negative dispersion bryk:jcp:2010 ; trachenko:rpp:2016 .

We conclude our investigation of the collective properties by plotting in Fig. 8 the density-density, density-energy, and energy-energy correlation functions at $T=0.10$ for two different wavevectors.

It is seen that these three different correlation functions have the same characteristics as discussed above, in agreement with the hydrodynamic predictions. That is, the standard DPD system includes the cross coupling between the longitudinal quantities hydrodynamically, at least qualitatively.

V Summary and Conclusions

In this paper the equilibrium relaxations of the standard dissipative particle dynamics model propopsed by Groot and Warren groot:jcp:1997 were investigated. First, the well-known results that the structure and dynamics at high temperatures ($T\geq 0.8$) resemble those of a gas were recaptured; this region in phase space is accordingly denoted the kinetic regime espanol:pre:1999 . At lower temperatures the viscosity increases with decreasing temperature and the Schmidt number approaches that of the model Lennard-Jones liquid. The DPD model features a single element Maxwellian shear modulus relaxation behavior for sufficiently small frequencies that depend on the temperature; the lower the temperature the smaller the frequencies that are required to observe Maxwellian behavior. Also, the time-temperature superposition principle is applicable in the low frequency regime.

For nonzero wavevectors, the hydrodynamic prediction for the transverse velocity autocorrelation function is tested using the Jung-Schmid and the Irving-Kirkwood definitions of the viscosity; the former includes the random and dissipative shear force contributions whereas the latter only includes the conservative and kinetic contributions. Using the Irving-Kirkwood viscosity the hydrodynamic predictions are in excellent agreement with simulations results for temperatures $0.4\leq T\leq 0.7$ and $0.53\leq k_{y}\leq 2.12$. Importantly, using the Jung-Schmid viscosity overestimates the relaxation, indicating that the transverse relaxation dynamics are independent of the dissipative and random shear forces. Also, for higher temperatures the agreement is less satisfactory, for a given wavevector, in accordance with the Bocquet-Chaix criterion.

A qualitative investigation into the longitudinal dynamics was also carried out. For the high temperature regime ( $T\geq 0.8$), the density longitudinal spectrum at low wavevectors is characterized by a single sharp Brillouin peak. This indicates that the longitudinal relaxation is athermal and dominated by propagating damped density waves. This mechanism is very different compared to a simple liquid, in which the thermal diffusion process dominates at low wavevector. In the low temperature range, the Rayleigh peak is more prominent; a fingerprint of the thermal diffusion process. Dispersion relations for the Landau-Placzek ratio shows that the thermal process intensity increases compared to the wave propagation process as the length scale decreases; this is true for all temperatures and wavevectors studied and also the case for the Lennard-Jones liquid, even though the Landau-Placzek ratio is larger here. For the supercritical fluid Lennard-Jones model there is a small increase in the speed-of-sound with respect to temperature bryk:jcp:2010 , however, for the DPD model this a constant with respect to temperature for $k<2$. Finally, the DPD model features the cross couplings predicted by the theory, at least, qualitatively.

In conclusion, the thermal fluctuations in the standard coarse grained DPD model by Groot and Warren groot:jcp:1997 preserves, at least qualitatively, the underlying mechanical processes predicted by classical hydrodynamic theory. Therefore, the model can be used to study fluctuating hydrodynamics as stated by Espańol and Warren espanol:jcp:2017 . However, we suggest to use low temperature settings where $T\leq 0.7$.

VI Acknowledgments

Innovation Fund Denmark supported this work as a part of the ROSE project (no. 5160-00009B). MLG acknowledges support from the US Fulbright Program. This work was supported by the VILLUM Investigator grant for the Matter project (JCD). Finally, we wish to thank the reviewers for their valuable comments.