Green-Kubo measurement of liquid-solid friction in finite-size systems

To characterize liquid-solid friction using molecular dynamics simulations, Bocquet and Barrat (BB) [Phys. Rev. E 49, 3079–3092 (1994)] proposed to use the plateau value of a Green-Kubo (GK) integral of the friction force. The BB method is delicate to apply in finite-size simulations, where the GK integral vanishes at long times. Here, we derive an expression for the GK integral in finite-size systems, based on a Langevin description of a coarse-grained system effectively involving a certain thickness of liquid close to the wall. Fitting this expression to GK integrals obtained from simulations of a liquid slab provides the friction coefficient and the effective thickness of the coarse-grained system. We show that the coarse-grained system for a Lennard-Jones fluid between flat and smooth solid surfaces is 2–3 molecules thick, which provides a criterion for measuring the friction coefficient independently of confinement. As compared to nonequilibrium simulations, the new approach is more accurate and removes some ambiguities of nonequilibrium measurements. Overall, we hope that this new method can be used to characterize efficiently liquid-solid friction in a variety of systems of interest, e.g., for nanofluidic applications.


I. INTRODUCTION
To describe flows in macroscopic channels, one can combine a continuum description of mass transport in the bulk (e.g., the Navier-Stokes equation) with a no-slip boundary condition (BC) for the fluid velocity at the wall. Such a macroscopic description is expected to fail to capture flows in nanofluidic systems, 1,2 which offer great promises of application in green energies [3][4][5] and water treatment. [6][7][8] Yet, it has been discussed that for water the continuum description of the bulk flow should remain valid down to typically 1 nm. 1 In contrast, both experiments and molecular simulations have shown that the no-slip BC can fail when liquid-solid friction is low enough. 9,10 In the presence of slip, the hydrodynamic BC is obtained by writing that the viscous shear stress in the liquid at the wall, η∂zv (with η being the shear viscosity, z being the direction normal to the interface, and v being the velocity parallel to the interface), is equal to the interfacial friction shear force per area τw, proportional to the slip velocity vs (the jump of parallel velocity at the interface), as initially discussed by Navier, 11,12 τw = λvs. (1) Equation (1) defines the (fluid) friction coefficient λ, a critical parameter controlling flows in nanofluidic systems. Molecular dynamics (MD) simulations, which explicitly describe the atomistic details of liquid-solid interfaces, can be instrumental in measuring the friction coefficient [13][14][15][16][17][18][19][20][21][22][23][24][25] and in exploring the underlying mechanisms. [26][27][28][29][30][31][32][33][34][35] Friction can be measured directly in nonequilibrium MD (NEMD) simulations of a flowing liquid.
However, the direct approach is time consuming and poses different issues 20,[36][37][38] related to, e.g., the choice of the thermostatting method, the necessary tests on the linear response of the system, or the definition of the effective position of the hydrodynamic boundary. Alternative approaches have been proposed based on equilibrium molecular dynamics (EMD). In particular, Bocquet and Barrat (BB) introduced a Green-Kubo (GK) expression where the friction coefficient is obtained as the plateau value at long times of a GK integral, 13 (2) where S, k B , T, and Fw denote the surface area, Boltzmann constant, absolute temperature, and instantaneous shear force exerted on the solid, respectively, and the angular brackets ⟨⟩ mean the ensemble average.
However, applying the BB formula in MD simulations poses a number of issues. 17,[21][22][23] In particular, simulations are limited to finite-size systems, where it has been discussed that GK integrals should vanish at infinite time instead of reaching a plateau corresponding to the response coefficient, 39,40 an issue sometimes referred to as "the plateau problem." 41,42 Different solutions have been suggested to the plateau problem. For instance, the maximum of the GK integral-or the value of the GK integral at the first zero of the autocorrelation function-has been used to estimate the friction coefficient. 13,[43][44][45][46][47] Alternative equilibrium [18][19][20]22,23,25,42 and nonequilibrium 14,15 methods have been developed. Recently, Español, de la Torre, and Duque-Zumajo 41 have proposed a general solution to the plateau problem within the context of Mori projection operator formulation, which could in principle be applied to the GK measurement of the liquid-solid friction coefficient. The fundamental statistical physics underlying the BB formula has also been discussed in two recent theoretical articles. 48,49 Nakano and Sasa 48 started from the Hamiltonian of the particle system and, by considering the local detailed balance condition, examined the applicability of existing expressions of the liquid-solid friction from a viewpoint of time-and length-scale separation, which can be evaluated from hydrodynamics behavior. On the other hand, Camargo et al. 49 considered the mechanical balance and local constitutive equation of a thin slab of layered fluid formed near the wall covering the liquid-solid interface (mentioned as a "pillbox" for a spherical interface) and derived the boundary condition. In contrast with these recent articles, here we propose a pragmatic approach to extend the BB formula to finite-size systems by introducing a simple expression for the memory kernel and fitting the full GK curve based on a coarse-grained description without dealing with the details of local mechanics. Based on this approach, we show the following two features: (1) an effective mass or thickness of liquid involved in the Langevin equation is obtained as a result and (2) the simple approach of taking the maximum of the GK integral leads to accurate measurements when the memory kernel time scale and the relaxation time scale as a function of the friction coefficient and effective mass are well separated. The two features could consequently correspond to the time-and length-scale separation mentioned in Nakano and Sasa, 48 and the effective coarse-grained system could also be related to the thin slab in Camargo et al. 49

II. THEORY
The BB formula can be derived based on a generalized Langevin equation (GLE), where M and U are the mass and velocity (relative to the wall) of a coarse-grained system of interest, and S, ξ, and R(t) denote the surface area, memory kernel, and random force, respectively. In their analytical derivation, Bocquet and Barrat 21 described the Brownian motion of a finite-mass wall in contact with a semi-infinite liquid, and M and U were the mass and velocity of the wall, respectively. In MD simulations, however, a liquid confined between two parallel immobile walls is usually considered. In that case, the coarsegrained system should effectively describe the hydrodynamic fluctuations of the liquid. We anticipate that the coarse-grained system can be described as an effective thin region of liquid close to the interface. However, we emphasize that the coarse-grained system encompasses the whole liquid dynamics and that the effective representation by a thin slab of fluctuating liquid should not be taken literally. In particular, viscous friction in the whole liquid is described implicitly through the coarse-graining process and we did not describe it explicitly in the following GLE description; this approach will be validated a posteriori by the simulations. We also emphasize that our model differs from previous work 18,19 where the equation of motion was applied to a liquid slab near the solid because we do not impose any constraint on the mass M of the coarse-grained system, which will be measured from the GK integral. A full expression for the GK integral can be derived from the GLE, Eq. (3), through the Laplace transform. The complete derivation can be found in the Appendix, and here we only summarize the key steps. Let the Laplace transform of a function f (t) be denoted by L(f ) ≡f (s). The Laplace transformCU (s) of the equilibrium autocorrelation function CU (t) of the velocity U writes On the other hand, the force Fw(t) exerted on the coarse-grained system of interest is given by Then, the autocorrelation function CF(t) of the force Fw satisfies the following relation: 50 Accordingly, the GK integral Λ(t) in Eq. (2) writes where dC U (t) dt | t=0 = 0 considering that CU (t) is an even function due to the stationarity condition. [51][52][53][54] Note that should be applied in the case the memory kernel ξ(t) is equal to the Dirac delta function, i.e., Eq.
dt has a discontinuous derivative at t = 0, see Refs. 21 and 54. The Laplace transform of Eq. (7) is Inserting Eq. (4) into Eq. (9), it follows that We now assume a simple Maxwell-type memory kernel 50 This memory kernel is originally designed to describe a viscoelastic response; here, we choose it for its simple functional form, although simulations will later show that tm is indeed compatible with viscoelastic relaxation. We also introduce the characteristic decay time of the GLE Note that through the coarse-graining process, hydrodynamic relaxation in the liquid-controlled by the kinematic viscosity νis enclosed implicitly in the coarse-grained mass M/S, and the decay time only depends explicitly on the liquid-solid friction coefficient.
Inserting the Laplace-transformed memory kernel in Eq. (10) and going back to the time domain, one can show that when t d > 4tm, the GK integral Λ(t) in Eq. (2) is given by where the 3 parameters λ 0 , t 1 , and t 2 can be obtained through 1 Equations (13)-(16) provide a generic framework that can be used to extract the friction coefficient from a fitting of the GK integral.
In particular, the friction coefficient can be obtained from the three fitting parameters, with u = t 2 /t 1 . One can also relate the maximum of the GK integral to λ and to the two time scales of the fit, From Eq. (18), it is clear that the maximum of the GK integral Λmax is not in principle equal to the friction coefficient λ and that Λmax and λ cannot be related without the knowledge of the characteristic times t 2 and t 1 , which can only be obtained through a full fit of the GK curve. However, when the time scales are well separated, i.e., in the limit of u → 0, Λmax becomes an increasingly good estimate of λ: the difference between Λmax and λ is on the order of 25% for u = 1/5 and falls, e.g., to 4% for u = 1/100 or to less than 1% for u = 1/1000. Note that in the limit where the time scales are separated, the two time scales t 1 and t 2 are given by t 1 ≃ t d and t 2 ≃ tm, respectively. Finally, in the limit of t d → ∞, so that the GK integral reaches a plateau equal to the friction coefficient, which corresponds to the original BB result.
In the following, we will test the general method presented here when the time scales of the memory kernel and of the GK integral decay are not well separated, using MD simulations.

III. SIMULATIONS
We considered a generic Lennard-Jones (LJ) liquid confined between parallel walls, see Fig. 1. We used fcc walls made of 8 atomic layers exposing a (001) face to the liquid; first neighbors in the solids were bound by a harmonic potential Φ h (r) = k/2 (r − req) 2 , with r being the interparticle distance, req = 0.277 nm, and k = 46.8 N/m. Interactions between liquid particles and between liquid and solid particles were modeled by a LJ pair potential, Φ LJ (r) = 4εij[(σij/r) 12 − (σij/r) 6 ], where i, j can be L for liquid particles or S for solid particles. The LJ interaction was truncated at a cut-off distance of rc = 3.5 σ and quadratic functions were added so that the potential and interaction force smoothly vanished at rc. 55 We used σ LL = 0.34 nm, ε LL = 121 K k B , σ LS = 0.345 nm, and ε LS = αε LL , where the "wetting coefficient" α was varied between 0.120 and 0.359. The corresponding contact angle of droplets with these parameters was shown in our previous work. 56 The atomic masses were m L = 39.95 u and m S = 195.1 u. We used periodic boundary conditions along the x and y lateral directions, with a box size Lx = Ly = 6.27 nm. The total system height (walls included) along the z direction was varied between 4 and 12 nm, with a corresponding number of liquid particles varying between ∼800 and 6400. We defined the liquid film height as H liq = z top wall − z bottom wall , with z top wall and z bottom wall being the positions where the liquid density vanishes, which can be accurately identified due to the steep rise of the density at the interface, see Fig. 1. Accordingly, H liq varied between 1 and 9 nm. The above definition of z wall was also used as the origin of δ eff described later in Eq. (20).
To compare the GK measurement and a reference NEMD measurement of the friction coefficient in the same system, we sheared the liquid in the x direction to extract the nonequilibrium response and we applied the GK formalism in the y direction. We checked on representative systems that independent EMD and NEMD simulations provided results consistent with the ones obtained using of Chemical Physics thermostatted at 100 K using a Langevin thermostat applied only in the z direction. We checked that the liquid temperature remained within 2 K of the set value. To set the pressure to 4 MPa, we first used the top wall as a piston, without shear; second, we applied shear, still using the top wall as a piston; and finally, we fixed the vertical wall position and continued shearing the system. We integrated the equations of motion using the velocity-Verlet algorithm, with a time step of 5 fs. The simulation time was 50 ns. Error bars in Figs. 2, 3, and 6 were obtained by separating the production run into five 10 ns chunks, measuring the quantities for each chunk, and estimating the statistical error within a 95% confidence interval from the five values.

IV. RESULTS AND DISCUSSION
We first computed the GK integrals Λ(t) for different wetting coefficients α, see Fig. 2. For low α (nonwetting surfaces), friction is low so that the decay time t d given by Eq. (12) is large and a plateau of the GK integral is observed at times intermediate between the memory kernel time tm and the decay time t d . However, as α increases, so does the friction coefficient. As a consequence, the GK integral goes to higher values, and the decay time decreases; for large α, t d becomes comparable with tm so that one cannot observe a plateau of the GK integral anymore. In that situation, the original BB formula of Chemical Physics ARTICLE scitation.org/journal/jcp cannot be directly applied and we used our finite-size extension of the formula, Eq. (13), to fit the GK integrals. Since at long times the error on the GK integrals becomes large, the numerical data were fitted for times t < 5 ps. The model indeed reproduces well the MD data at short and intermediate times, see Fig. 2. From the fit of the GK integrals, one can extract the values of the friction coefficient λ fit , of the coarse-grained mass per unit surface M/S, and of the memory kernel relaxation time tm. The evolution of these parameters with the wetting coefficient α is shown in Fig. 3. While the friction coefficient strongly depends on the wetting coefficient α, the memory kernel time tm remains constant within the error bars, at a value of ∼150 fs (except at the highest α where it is larger by ≲10%). This value is typical of viscoelastic relaxation in liquids. The effective mass of the coarse-grained system in the Langevin description increases slightly beyond the error bars with α, from ∼1 to 1.5 mg/m 2 . The fact that the effective mass does not depend much on α suggests that the description of the decay time as a function only of the liquid-solid friction coefficient and of a coarse-grained mass (enclosing effectively the liquid hydrodynamic modes), Eq. (12), is adequate in all systems considered. Note however that for the largest values of α, corresponding to the largest friction coefficients and smallest decay times t d , a slower decay of the GK integral seems to appear at long times, which is not captured by Eq. (13). This slower time scale could be related to momentum diffusion through the liquid, which by construction the coarsegrained description cannot capture. We will come back to that point when discussing the dependency of the GK integral with the system height.
We then computed the effective thickness δ eff of the coarsegrained system corresponding to the effective mass M/S, by integrating the liquid density profiles ρ liq (z) near the wall, In Fig. 3, one can observe that the increase of δ eff with α is similar to that of M/S. Interestingly, δ eff is on the order of 2-3 molecular sizes, i.e., the coarse-grained system of interest in the Langevin equation involves only a thin region of liquid close to the interface. Accordingly, the GK integrals should not depend on the system height for heights larger than δ eff .
We tested the influence of system height on the GK integral for different wetting coefficients α. Figure 4 shows that, as expected from the small effective thickness of liquid whose fluctuations control the fluctuations of friction force, the GK integrals obtained for different liquid heights H liq above 2.17 nm overlap so that their fit provides the same values for λ fit , M/S, and tm. On the other hand, the GK integrals obtained for the smallest liquid height H liq = 1.15 nm, which is comparable to δ eff in Fig. 3, are significantly different from those for H liq ≥ 2.17 nm. This indicates that a certain liquid thickness is needed to apply the present method.
More quantitatively, Fig. 5 shows the dependence of the friction coefficient, memory time, effective mass, and corresponding thickness on the liquid height H liq for systems with α = 0.240. As suggested by the GK curves in Fig. 4, neither λ nor M/S depend on H liq above a certain value so that the decay time t d = M/(λS) is also independent of the liquid height. The independence of t d on the system size seems to suggest that the GK integral Λ(t) will always vanish at infinite time regardless of the system size, and in particular at the thermodynamic limit. However, we have already mentioned that the coarse-grained description, while it worked very well to describe the GK integral at short and intermediate times, failed to describe a slower decay of the GK integral at long times. One can imagine two other characteristic time scales related to the momentum diffusion, i.e., t visc (δ eff ) ∼ δ 2 eff /ν and t visc (H liq ) ∼ H 2 liq /ν, where ν is the kinematic viscosity of the liquid. Using the bulk value ν ≈ 1.1 × 10 −7 m 2 /s, the former is roughly estimated about 9 ps with δ eff ∼ 1 nm, and the latter is from 12 ps to 670 ps with H liq for the present systems. These longer time scales should largely affect the asymptotic behavior. In particular, the second time scale diverges for infinite system height.
Finally, we compare the presented method to estimate the friction coefficient, λ fit , with the reference NEMD measurement, and with the maximum of the GK integral, Λmax. The NEMD friction coefficient was computed using Eq. (1), where we tested two definitions for the slip velocity: first, we measured the difference between the wall velocity and the liquid velocity at the position of the first adsorption layer z 1 , vs = |v liquid (z 1 ) − v wall |; second, we computed the difference between the wall velocity and the extrapolated bulk liquid velocity profile at the position of the first adsorption layer, vs = |v extrapolated (z 1 ) − v wall |. The NEMD friction coefficients obtained from these two definitions are denoted as λ NE,liq and λ NE,ext , respectively. Figure 6 compares the different measurements for varying the wetting coefficient. First, we note that the difference between λ fit and Λmax is less than 5%. Since in the systems considered, the time scale ratio u = t 2 /t 1 was on the order of 1/100 or less, the relative error between Λmax and λ fit is consistent with the estimate provided in Sec. II. In practice, this means that, at least for the systems considered here, simply measuring the maximum of the GK integral provides a good estimate of the friction coefficient. We then turn to the NEMD measurements. One can notice that the error bars are much larger for the NEMD points than for the EMD points, even though the data were obtained for the same simulation time. This suggests that EMD provides a more efficient route to measure liquid-solid friction than NEMD. In addition, two different approaches to measure the slip velocity, i.e., based on the fluid velocity in the first adsorption layer, λ NE,liq or based on the extrapolated bulk velocity profile at the same position, λ NE,ext , result in different values for the friction coefficient. The NEMD measurement will also depend on the position where the slip velocity is measured. While using the position of the first adsorption layer is common practice, this approach has no firm fundamental basis.
Overall, the NEMD estimate is therefore less accurate for a given simulation time and sensitive to the procedure used to determine the slip velocity. In contrast, fitting the GK integral provides an accurate and unambiguous value. Although the EMD and NEMD seem to deviate at large wetting coefficients α, the difference is within the error bars. Additionally, we would like to emphasize that for the largest values of α, the friction coefficient is very large and corresponds to very small slip lengths, b ∼ 1 nm. For such values of b, any uncertainty on the hydrodynamic wall position will result in a large error on the friction coefficient so that we suggest that EMD should here be thought as the reference measurement method and that the apparent discrepancy between EMD and NEMD highlights the difficulties of getting a proper NEMD estimate in the studied systems.

V. CONCLUSION
The BB formula, identifying the liquid-solid friction coefficient with the plateau value of a GK integral, poses some problems in finite-size simulations, where the GK integral vanishes at long times. This formula can be obtained from a Langevin description of a diffusing wall in contact with a semi-infinite liquid. Here, we derived the analytical expression of the GK integral for a finite-size liquid slab confined between immobile walls by applying the Langevin description to a coarse-grained system involving effectively a fraction of the confined liquid, assuming a simple Maxwell-type memory kernel relating the slip velocity and velocity of the coarse-grained system. Using this generic analytical framework, we showed that a common procedure of taking the maximum of the GK integral as the friction coefficient was reasonably accurate when the memory time and the decay time of the GK integral were well separated (e.g., a 4% error for a factor of 100 between the two times).
By fitting MD results to the derived expression, we could extract the friction coefficient, the memory time, the mass of the coarsegrained system, and the corresponding effective thickness of liquid whose fluctuations control the fluctuations of the friction force. We varied the wetting properties by tuning the liquid-solid interaction potential. The friction coefficient was strongly affected by wetting, while the memory time remained approximately constant, at a value consistent with viscoelastic relaxation. The effective thickness δ eff of fluctuating liquid did not change much and was on the order of 2-3 molecular sizes. Accordingly, we checked that the measured friction coefficient was independent of the liquid film height H liq when it remained large as compared to δ eff . For a smaller system with H liq ∼ δ eff , the GK integrals changed. The estimate of δ eff from the fit of the GK integral therefore provides a self-consistent criterion on of Chemical Physics ARTICLE scitation.org/journal/jcp the minimal system height that needs to be used to characterize interfacial friction without confinement effects. Finally, we compared EMD and NEMD measurements. For a given simulation time, the approach presented here provides a lower statistical error; additionally, in NEMD an ambiguity on the measured friction coefficient results from the difficulty of defining the slip velocity, but in the proposed method the slip velocity does not need to be defined. Overall, we suggest that the proposed approach is simple and accurate, and provides an efficient path to characterize liquid-solid friction in MD simulations.

APPENDIX: COMPLETE DERIVATION
Here is the derivation of Bocquet-Barrat formula We start from the generalized Langevin equation (GLE) 51 where M and U are the mass and velocity of the coarse-grained system of interest, respectively, while S, λ ′ , and R(t) denote the surface area, retarded effect of the friction, and random force, respectively. The macroscopic friction coefficient λ is related to λ ′ by 45 By substituting λ ′ (t) = λξ(t) with a memory kernel ξ, which satisfies it follows for Eq. (A2) that Note that in the case the wall is mobile as in the theoretical derivation of Bocquet and Barrat, 21 the mass and velocity are replaced by wall mass and velocity Mw and Uw, respectively, and the slip velocity vs and Uw are related by A Green-Kubo relation is derived from the GLE, Eq. (A5), through the Laplace transform. By multiplying U(0) to both sides of Eq. (A5) and by taking the ensemble average, it follows that Denoting CU (t) ≡ ⟨U(0)U(t)⟩ the equilibrium autocorrelation function of U and assuming ⟨U(0)R(t)⟩ = 0 for t ≥ 0, Eq. (A7) is rewritten as where the lower limit of the integral is changed from −∞ to 0 so that CU (t) satisfies the stationarity condition, 51,52 i.e., CU (t) is an even function. In addition, the decay time scale t d given by On the other hand, the force Fw(t) exerted on the coarse-grained system of interest is given by Then, the autocorrelation function CF(t) = ⟨Fw(0)Fw(t)⟩ of the force Fw, satisfies the following relation: 50 By using this, we examine the function form of the integrated force autocorrelation function Λ(t) given by where dC U (t) dt | t=0 = 0 considering that CU (t) is an even function due to the stationarity condition. [51][52][53][54] Note that of Chemical Physics ARTICLE scitation.org/journal/jcp should be applied in the case the memory kernel ξ(t) is equal to the Dirac delta function, i.e., Eq. (3)