Molecular hydrodynamics from memory kernels

The memory kernel for a tagged particle in a fluid, computed from molecular dynamics simulations, decays algebraically as $t^{-3/2}$. We show how the hydrodynamic Basset-Boussinesq force naturally emerges from this long-time tail and generalize the concept of hydrodynamic added mass. This mass term is negative in the present case of a molecular solute, at odds with incompressible hydrodynamics predictions. We finally discuss the various contributions to the friction, the associated time scales and the cross-over between the molecular and hydrodynamic regimes upon increasing the solute radius.

The Brownian motion of a particle in a fluid finds its origin in the fluctuating force exerted by the solvent molecules on the solute. It has long been known that the canonical description of this random force by a Gaussian Markov process is only valid in limiting cases. Even in the limit where the solute is much heavier than the solvent particles, for which multiple time-scale analysis allows to recover the Smoluchowski equation for diffusion [1], non-Markovian effects are expected when the mass density ratio is close to unity [2] -a situation which is rather the rule than the exception e.g. in colloidal suspensions. These non-Markovian effects arise because of momentum conservation, leading to slow hydrodynamic modes that manifest themselves as long-time tails in the velocity autocorrelation function (VACF) [3][4][5][6]. Recent experiments have demonstrated that the force exerted by the bath includes a deterministic component [7], well described for large colloidal spheres by the Basset-Boussinesq (BB) hydrodynamic force [8,9]: where R is the sphere radius, η the solvent viscosity and ρ 0 its mass density. The first term is the usual Stokes friction. The other two account for the inertia of the displaced fluid and involve a finite added mass m BB 0 = 2 3 πR 3 ρ 0 and a viscosity-dependent retarded component describing the transient effects of momentum diffusion in the solvent. While continuous descriptions of steady-state flows appear to hold down to the nanoscale [10][11][12], possibly at the price of adapting the hydrodynamic radius R or the boundary conditions [13], their validity for the transient regimes should be questioned. The implicit assumption of a separation of time scales between the solvent and solute dynamics, which holds a priori for colloidal parti-cles [14], is expected to break down with smaller solutes such as nanoparticles or biomolecules.
Here we address the fundamental questions that arise when approaching the regime of molecular solutes by computing directly from Molecular Dynamics (MD) simulations the memory kernel and the random noise of the Generalized Langevin Equation (GLE). A novel algorithm based on the Mori-Zwanzig formalism with high numerical stability allows us to explore long time scales for the first time. We consider the extreme case of a tagged particle (identical masses and sizes) in a pure supercritical fluid.
By examining the long-time behaviour of the memory kernel, we demonstrate the generality of the functional form of Equation 1 beyond pure hydrodynamic descriptions and discuss its interpretation as the time-dependent force exerted by the solvent on the solute at thermal equilibrium. Importantly, we show how to define and compute a mass from the memory kernel itself. This generalisation from the microscopic dynamics correctly describes the numerical results for the VACF almost down to the ballistic time scale and provides insights into the emergence of the hydrodynamic behaviour for larger solutes, bridging the gap between the solvent and colloidal time scales.
In the Zwanzig-Mori formalism [15][16][17], the velocity v(t) of a tagged particle of mass m in a fluid follows the generalized Langevin equation where K(u) is the memory kernel and R(t) the so-called random force, that are obtained from the true force F acting on the tagged particle using the projection operator technique and defined as with k B the Boltzmann constant and T the temperature, and R(t) = e i(1−P)Lt F. In these equations, iL is the Liouvillian operator corresponding to the unperturbed dynamics and P is the Mori projection operator along the velocity v, acting on an observable A as PA = vA v 2 v. Throughout the paper, · denotes the canonical equilibrium average at temperature T . The force F is propagated using the orthogonal dynamics e i(1−P)Lu instead of the normal dynamics to obtain the Zwanzig-Mori memory kernel.
The auto-correlation function of this projected force, or noise, F e i(1−P)Lu F differs significantly from the autocorrelation function of the force F F(u) . In particular, for a periodic system the latter integrates to zero whereas the former integrates to the friction ξ. This property of the projected force is a form of the Einstein relation since the friction is related to the diffusion constant D of the tagged particle by D = kB T ξ . However, extracting the projected force correlation function, or the kernel, from MD simulations is a difficult task.
We have recently introduced two practical schemes to compute such properties for generic observables from MD trajectories [18]. These algorithms are only accurate to first order in the MD timestep δt -thus preventing their use to investigate the long time behaviour. Here we employ a novel algorithm [19], which provides second order accuracy at virtually no additional computational cost, to study the memory kernel for diffusion in a Lennard-Jones (LJ) fluid. We consider a system of 10 4 LJ particles at a reduced density ρ * = ρσ 3 = 0.5 and reduced temperature T * = k B T /ǫ = 1.5 with σ and ǫ the LJ diameter and energy, respectively, i.e. at the critical density and slightly above the critical temperature. Newton's equations of motion are solved using the velocity Verlet algorithm and cubic periodic boundary conditions. Interactions are computed using a cut-off radius r c = 3σ. The system is first equilibrated at the target temperature during 230.41 t * by performing MD with a timestep of 9.2×10 −4 t * , in the NVT ensemble using Langevin thermostat with a time constant of 0.92 t * . All properties are then determined from a 230.41 t * trajectory with a timestep of 4.6×10 −4 t * in the NVE ensemble generated with the DLPOLY [20] simulation package and block averages were taken over trajectory segments of one tenth of the total trajectory.
The novel second order algorithm presents remarkable long time stability and allows to investigate time scales much beyond ∼ t * . This is demonstrated in the inset of the figure 1, which displays the running time integral of the noise auto-correlation function (NACF). From the plateau of the NACF (Fig. 1), we obtain ξ ∼ 4.5 ± 0.1 LJ units, in excellent agreement with the Einstein relation (k B T /D ∼ 4.4 LJ units). In contrast, the running time integral of the unprojected force auto-correlation function (FACF) tends to zero as expected.
where ν = η ρm is the kinematic viscosity, with η the fluid viscosity and m the particle mass. The diffusion constant D is often omitted in this long time tail, however it is necessary to reproduce our numerical result as can be seen from figure 1. This term is due to the diffusion of the particle simultaneously with the momentum transfer in the fluid [21,22]. The FACF is the second-order derivative of the VACF and should decay in the same limit as: This is indeed the case as shown in figure 1. In contrast, the NACF, which is nothing but the memory kernel K, decays much more slowly than the FACF, following the same t − 3 2 scaling as the VACF. In fact, such a scaling is not unexpected: Corngold indeed showed from the relation between the Laplace transforms of Z(t) and K(t) that under rather mild conditions for the VACF, the memory kernel defined by Eq. 2 should decay as [23] K(t) ∼ − ξ 2 kB T Z(t) leading to from the asymptotic behavior of Z(t). As can be seen in figure 1, this prediction is indeed satisfied by the memory kernel determined from MD. This scaling is also consistent with the low frequency limit of the hydrodynamic memory kernel corresponding to Eq. 1 (see below). It has then been observed experimentally for colloidal particles where this limit applies [7]. Our results confirm for the first time that this scaling also holds for the diffusion of microscopic particles. From the decay of the memory kernel at long times, the development of the Laplace transform of the friction kernel isK with: where we have introduced a mass defined by: under the assumption that K(t) + 1 2 απ −1/2 t −3/2 decreases to zero faster than t −2 , and where an integration by parts was used for the second equality. Note that while the speed of convergence depends on higher order terms in the expansion Eq. 7, the value of m 0 defined by Eq. 9 does not. Figure 2 shows the running integral associated with the definition of this mass. The observed plateau demonstrates the convergence of the integral and thus validates the above assumption in the present case. Eq. 9 therefore provides the first definition of the mass term from the microscopic dynamics. Surprisingly, this mass term is negative, with a value of m 0 ∼ −0.18 m -in contradiction with the incompressible hydrodynamic prediction for the added mass m BB 0 = 2 3 πR 3 ρ 0 [22,24]. This observation can be interpreted as follows, by analyzing the various contributions to the kernel K. At short times, K is dominated by short-range collisions between the solute and the solvent and can be approximated by an exponential decay K(t) = ξE τ0 e −t/τ0 , with ξ E the Enskog friction [4] and τ 0 the characteristic time for the decay of the FACF (τ 0 ∼ 0.05 t * in the present case). The subscript 0 indicates that this time corresponds to the collisions between the solvent molecules and the solute, rather than to a time scale associated with the decay of the solute VACF. This collisional component of the kernel contributes to the mass defined in Eq. 9 as a negative term Computing the Enskog friction [25] for a solute of size σ, we get ξ E ≈ 5.8. This value is consistent with the maximum of the time-dependent friction in figure 1 (see below) and results for the mass to a contribution m E 0 ≈ −0.29 m. Other mechanisms come into play on times scales longer than τ 0 . Indeed, momentum transfer from the solute to the solvent includes a transient regime giving rise to a positive contribution to the mass term (here over a time τ m ∼ 0.5 t * , as can be seen in Fig. 2) and eventually becomes diffusive, leading to the retarded force and to a decrease in the friction (see Fig. 1): The solvent backflow tends to drag the solute in the direction of its initial velocity, i.e. contributes negatively to the friction. Assuming that this component of the mass term is well described by the hydrodynamic result despite the molecular size of the solute, we obtain a total mass m 0 = m E 0 + m BB 0 ≈ −0.16 m, which is in good agreement with the MD result considering the strong assumptions involved (validity of the Enskog result at high packing fraction and hydrodynamic model of the mass), and conforts our interpretation of the two competing contributions to the mass term.
We now consider the ensemble-average velocityv obtained over an ensemble of identical systems initially in equilibrium and put out of equilibrium at time t = 0 by a time-dependent applied force f ǫ (t), identical to all replicas of the system. We show in Supplementary Material that the evolution of the ensemble average velocitȳ v is given by the same kernel as the GLE for the microscopic velocity with the random force replaced by the applied force [19]. For slowly varying forces, the ensembleaverage velocity also varies slowly and we can consider the s → 0 limit inK(s). The first three terms of Eq. 7 correspond to an evolution ofv, according to: for a system put out of equilibrium from t = 0 by a slowly-varying infinitesimal force f ǫ (t). There is no hypothesis of separation of time scales between slow and fast degrees of freedom of the system in this equation, its meaning is that of a slowly varying response to a slow perturbation [26]. This evolution provides a generalisation of the Basset-Boussinesq Eq. 1 to arbitrary solutes satisfying only the above generic assumptions on the long-time behaviour of the corresponding memory kernel. Following the method of Chow and Hermans for the VACF of a particle subject to a BB force [9], we express analytically the VACF of the solute subject to the force Eq. 10 and compare it to the simulation results in Figure 3, both in logarithmic and linear scales. The agreement is excellent down to relatively short times (less than 0.5t * ), without any adjustable parameter. This further demonstrates the relevance of the above definition of the mass term from the memory kernel (i.e. not from the solute geometry and hydrodynamic properties of the solvent). Note, however, that truncating the memory kernel to the first three terms of the low frequency expansion Eq. 7 leads to some limitations for the description of the short-time behaviour, such as an incorrect initial value of the VACF, namely k B T /(m + m 0 ) instead of k B T /m [9].
A negative contribution to the mass term can also be derived from the hydrodynamics of compressible fluids, involving the time R/c it takes for sound waves to propagate over the particle radius -confirming the role of retardation effects in this negative contribution. However, introducing compressibility in continuum hydrodynamics [27,28] does not improve the prediction for the VACF, even with an effective hydrodynamic radius adjusted to reproduce the calculated friction (see [19]), because it does not capture molecular scale effects.
FIG. 4: Schematic evolution of the friction (a) and mass (b) with increasing solute radius R for a spherical density-matched solute. The curves correspond to increasing radius from black to red to blue.
Finally, let us consider the implications of the present work for larger solutes. We consider here spherical solutes with a density equal to that of the solvent, which is the most common experimental situation of density-matched colloidal suspensions, with a mass M = 4πR 3 3 ρ 0 . The following discussion is illustrated in Figure 4. For a large particle (R ≫ σ/2 and M ≫ m), the Enskog friction ξ E ∝ R 2 and the negative Enskog contribution −ξ E τ 0 to the mass (where τ 0 only weakly depends on R) dominate at short times. The Stokes friction ξ S ∝ R and the BB hydrodynamic mass m BB 0 ∝ R 3 are recovered over a time τ m ∝ R 2 /ν. While in the present case this analysis neglects molecular features, τ m provides the correct order of magnitude for the time over which the integral defining the mass converges (see Fig. 2).
The long time-tail of the memory kernel bridges microscopic dynamics with continuum hydrodynamics as it gives rise to a force entering in the evolution equation of the tagged velocity similar to the BB hydrodynamic force. The memory kernel further allows for the first microscopic definition of the mass present in this evolution equation. This mass is found to be negative for a solute identical to solvent particles and is related to the retardation of the friction force. Extracting the mass term directly from MD simulations paves the way to the study of isotopic effects. It can also be used to quantify in a well-defined way the number of molecules brought along ions during transport or to interpret the peculiar behaviour of the friction on alcanes as a function of chain length [29,30]. In particular, it provides a microscopic route to model acoustophoresis [31,32] or electro-osmotic effects [33,34]. Finally, the novel algorithm introduced here could be used to compute projected correlation functions of other observables and to investigate the properties of the corresponding GLE.
The authors are endebted to Jean-Pierre Hansen and Lydéric Bocquet for very fruitful discussions and for critical reading of the manuscript.