Theoretical framework for the atomistic modeling of frequency-dependent liquid-solid friction

Nanofluidics shows great promise for energy conversion and desalination applications. The performance of nanofluidic devices is controlled by liquid-solid friction, quantified by the Navier friction coefficient (FC). Despite decades of research, there is no well-established generic framework to determine the frequency dependent Navier FC from atomistic simulations. Here, we have derived analytical expressions to connect the Navier FC to the random force autocorrelation on the confining wall, from the observation that the random force autocorrelation can be related to the hydrodynamic boundary condition, where the Navier FC appears. The analytical framework is generic in the sense that it explicitly includes the system size dependence and also the frequency dependence of the FC, which enabled us to address (i) the long-standing plateau issue in the evaluation of the FC and (ii) the non-Markovian behavior of liquid-solid friction of a Lennard-Jones liquid and of water on various walls and at various temperatures, including the supercooled regime. This new framework opens the way to explore the frequency dependent FC for a wide range of complex liquids.

Introduction-Nanofluidics is the discipline that describes fluid motion in nano-confinement, whose unique behavior in the mass and ionic transport should be a key ingredient in future technologies for fluid filtration and energy harvesting [1][2][3][4][5]. The recent advent of new materials and fabrication techniques to create fine fluid conduits [6] has even increased the importance to explore nanofluidic transport. In nanofluidic systems, surface effects play a critical role because of the large surface-tovolume ratio. In particular, liquid-solid slip can boost the performance of nanofluidic devices [7][8][9][10].
Liquid-solid slip was first foreseen by Navier [11], who proposed that the slip velocity u slip is proportional to the shear stress on the wall τ as τ = λ u slip (1) with λ the Navier friction coefficient (FC). The Navier boundary condition, Eq. (1), has been tested by many authors in the past decades as reviewed in the articles [2,[12][13][14]. From the theoretical side, the pioneering work by Bocquet and Barrat [15,16] showed that λ can be related to the equilibrium fluctuations of the friction force through a Green-Kubo (GK) formula: where δF is the random friction force on the wall at equilibrium, S the wall surface area and T the temperature, with k B the Boltzmann constant. Later, the system size dependence of the formula, called the plateau problem, was pointed out [17] and alternative ways to estimate the FC have been proposed [18][19][20][21].
Recently, two theoretical approaches have been proposed to challenge the plateau problem fundamentally. Español and coworkers [22,23] developed a new theory of non-equilibrium statistical mechanics, which led to a corrected form of the GK formula under the assumption that the system is Markovian. In another recent work, Nakano and Sasa [24,25] introduced explicit assumptions on the scale separation between the microscopic motion of molecules and the macroscopic motion of fluid and proposed a new way to estimate λ based on linearized fluctuating hydrodynamics. Both works from the two groups involved elaborate mathematical manipulations, and only reported the pure viscous (Markovian) behavior of the Navier FC, for a Lennard-Jones (LJ) liquid on a simple model wall. However, non-Markovian behavior of the FC was recently reported for a LJ liquid on a fcc lattice [26] and it is plausible that more complex liquids such as water also show such behavior, in analogy with their bulk transport properties [27][28][29][30][31][32].
In this Letter, we develop a theory to relate the Navier FC and the random force autocorrelation on the wall, by employing rather classical tools such as Stokes and Langevin equations. This theoretical framework offers some new perspectives on long-standing debates related to the GK modeling of liquid-solid friction, together with a simple and fast method to fully characterize the frequency-dependent Navier FC. We then apply this method to explore the frictional behavior of a simple LJ liquid and of water at various temperatures -including the supercooled regime.
Theory-Let us consider the system shown in Fig. 1(a), where a liquid is confined between two solid walls under no external field. When the bottom wall is let to move freely in a wall-tangential direction x, its motion can be described by a Langevin equation [33]: where M , S and U are the mass, the surface area and the x-direction velocity of the bottom wall, respectively; ξ is the friction kernel and δF bot is the random force that originates from the direct interaction between the solid and liquid particles. Assuming the equipartition of energy, Eq. (3) leads to the fluctuation-dissipation theorem: The motion of the liquid in response to the bottom wall motion (the top wall is fixed) can be described by Stokes equation for a wide frequency range [26] [34]: with the Navier boundary condition defined on the bottom and top hydrodynamic boundaries where u, t, ρ, η and λ denote the liquid velocity in the x direction, the time, the bulk liquid density, the bulk liquid viscosity, and the Navier FC, respectively. Note that λ is frequency dependent and of non-Markovian nature: λ has the dimension [Pa/m] instead of [Pa · s/m] adequate for the frequency-independent FC. Because the first term on the RHS of Eq. (3) would also be written ] dt ′ in terms of the slip velocity on the wall, the friction kernel ξ is given from the solution of Eqs. (5) and (6) (for the complete derivation, see the supplemental material, SM). Combined with Eq. (4), the expression for the force autocorrelation function is obtained: where the tilde indicates that the variables are Fourier-Laplace transformed and ζ denotes iρω/η, with ω the angular frequency. Considering that the nature of the random force is independent of the wall velocity by construction, Eq. (7) holds even when the bottom wall is fixed. Unlike the GK formula by Bocquet and Barrat [33], this equation explicitly includes the system size dependence in the relation between the random force autocorrelation and the Navier FC. It is also more general as it provides the viscoelastic behavior of the friction coefficient. From Eq. (7), the asymptotic behaviors of the random force autocorrelation are lim ω→0, h:finiteC δF bot where λ 0 is the zero-frequency component of the Navier FC and b is the slip length defined as η/λ 0 . From Eqns. (8) and (9), one can derive several important properties of the random force, whose evidence will be shown later in the Results section. First, the GK integral of the random force (Eq. 8) is not zero but has a finite value that depends on the system height h, which tells that the integral has a plateau. Note that the Bocquet-Barrat formula may be recovered by taking the plug-flow limit b ≫ h. In this limit only, the GK integral no longer depends on the system height h. Second, in the thermodynamic limit where the system height is infinite, the GK integral of the random force goes to zero: this is true regardless of the order in which the limits are taken, lim h→∞ lim ω→0C δF bot = lim ω→0 lim h→∞C δF bot = 0. Finally, when the frequency ω is high enough so that the penetration length is much smaller than the magnitude of the complex slip length ( η/ρω ≪ |η/λ|) as well as than the system height ( η/ρω ≪ h), the random force autocorrelation coincides with the Navier FC: for small t satisfying ηt/ρ ≪ min{|η/λ|, h}. All these results, which are in contrast to the common view that λ 0 might be obtained as lim ω→0 lim h→∞C δF bot [33,35], reflect our explicit consideration of the system height and its hydrodynamic influence on the fluctuations of the friction force. Interestingly, one can find the Navier FC for the whole frequency range from the measured random force autocorrelation C δF bot by solving Eq. (7) forλ, which is a quadratic equation of it. The physically correct solution out of the two can be chosen so that it satisfies the following relationship: where C δF total is the autocorrelation function of the random force summed over both the top and bottom walls (for the complete procedure, see the SM). Simulation-To validate the ideas in the Theory section, we performed equilibrium molecular dynamics (MD) simulations for two kinds of liquid on different walls.
The first system was a LJ liquid confined between two fcc crystal walls (Fig. 1b). The liquid consisted of 6400 molecules unless otherwise mentioned. The quantities for this system are shown in LJ reduced units based on the liquid molecular mass m f (= 6.634 × 10 −26 kg) and two parameters σ ff (= 3.40Å) and ε ff (= 1.67 × 10 −21 J) for the LJ potential Φ LJ (r ij ) = 4ε (σ/r ij ) 12 − (σ/r ij ) 6 between the liquid molecules, with r ij the distance between particle i and j with a cut-off distance of 3.50σ [36]. Each solid wall consisted of 8 layers of atoms in the (001) plane of a fcc crystal with a lattice constant of 1.15σ ff . The system temperature was controlled at 0.827ε ff /k B by a Langevin thermostat set on the second outermost layer of the walls, and the pressure was set to 0.094ε ff /σ 3 ff by a preliminary piston equilibration (see Ref. 26 for technical details). For the liquid-solid interaction, the LJ potential was adopted as well. To see the effect of the wettability, which is known to have impact on the friction [2], three different ε fw were used, ε fw = ε 0 fw = 0.155ε ff , ε fw = 2ε 0 fw , and ε fw = 3ε 0 fw , while σ fw = 1.01σ ff was kept constant. The corresponding contact angles of a sessile LJ droplet on the three walls were 136, 79 and 0 degrees, respectively [37].
The other two systems were water confined between either graphene walls or fcc crystal walls (Figs. 1c and 1d, respectively). In this case the fluid was constituted by 4096 TIP4P/2005 water molecules [38]. For water enclosed between generic fcc walls, such walls were constituted by three atomic layers of a fcc crystal exposing the (001) face with a lattice constant of 5.356Å, and with liquid-solid interaction parameters taken for hydrophobic walls from Ref. 39. For water enclosed between graphene walls, the liquid-solid interaction parameters were taken from Ref. 40. For both systems, the temperature was controlled by applying a Nosé-Hoover thermostat to liquid atoms, and the pressure was set to 1 atm through a preliminary piston equilibration (see Ref. 41

for technical details).
Results-First, we discuss the convergence of the GK integral of the random force autocorrelation. Figure 2 shows the normalized finite-time GK integral, Λ(t) := t 0 C δF bot (t ′ ) dt ′ /Sk B T , as a function of the upper time limit of the integration, for the LJ liquid confined by walls with ε fw = 2ε 0 fw . To obtain Λ(t) up to t = 1000 ps, we produced the simulation data typically for 300 ns. One can see that Λ(t) has a system-height-dependent plateau for t → ∞, whose value, the GK integral, de-creases by increasing the system height. The figure also illustrates that the whole Λ(t) profile is well reproduced with C δF bot calculated from the RHS of Eq. (7), and the plateau values coincide with the RHS of Eq. (8). In this evaluation of Eqs. (7) and (8), we substituted λ(t) by the Maxwell-type model λ 0 exp(−t/t λ )/t λ with the parameters λ 0 = 0.1492 √ m f ε ff /σ 3 ff and t λ = 0.077σ ff m f /ε ff taken from the results of non-equilibrium simulations [26], whose simulation system and conditions were identical to the present study [42]. Here, the same λ was used regardless of the system height: this shows that there is no system size dependence in the estimation of λ from the MD simulation data by the present theory. Sometimes in the literature [16,43], λ 0 is estimated as max Λ(t), which gives about 0.14 √ m f ε ff /σ 3 ff and slightly underestimates λ 0 . This under-estimation was also shown in Ref. 21. In summary, the discussion here provides a new perspective to the long-standing plateau issue for the evaluation of the GK integral: there is a plateau in the limit of the integral for finite-sized systems and this plateau value has a hydrodynamic meaning. Now we determine the complete nature of the Navier FC for the three liquid-solid interfaces shown in Fig. 1. As described in the Theory section, the Navier FC λ can be obtained by solving Eq. (7) for λ once the bulk liquid properties and the random force autocorrelation are measured. To estimate the hydrodynamic system height h, for the LJ liquid system we employed the values from Ref. 26 and for the water systems we adopted the separation between the Gibbs dividing surfaces [44] on the two confining walls. Figure 3 shows the Navier FC λ(t) and the normalized random force autocorrelation C δF bot (t)/Sk B T for the LJ liquid and for water. The equivalence between λ and C δF bot /Sk B T , Eq. (10), holds almost everywhere, although the condition ηt/ρ ≪ min{|η/λ|, h} is not strictly satisfied in the tail region for supercooled water. Figure 3(a) shows λ(t) for the LJ liquid system, together with a Maxwell-type viscoelastic model λ M (t) = λ 0 exp(−t/t λ )/t λ [26] for comparison. The Maxwell-type model describes well the time-dependent behavior of the FC for this simple liquid, which was expected from Fig. 2 showing good reproduction of the GK integral applying λ(t) = λ M (t) in Eq. (7). A slight distortion in the λ M profile comes from the non-vanishing time derivative at t = 0. Note that dλ/ dt| t=0 = 0 follows from Eq. (10) because the autocorrelation function is an even function in a stationary system. Figures 3(b) and (c) show λ(t) for water confined either by graphene walls or by fcc walls at three different temperatures: 235 K, 268 K and 360 K. The profiles of the FC on both walls at 360 K (i.e. for liquid water above its melting point) look similar to those of the LJ liquid. However, for supercooled water, i.e. metastable liquid water below its melting point, λ(t) cannot be described by the Maxwell-type model. It is known for supercooled liquids that the density relaxation is a two-step process with two characteristic decay times [45,46]: here one can see that the FC also decays with two characteristic times.
Conclusions-We have derived analytical expressions to connect the equilibrium fluctuations of the random force on the wall and the Navier friction coefficient (FC). The expressions are generic in the sense that they explicitly include the system size dependence and also the FC can be frequency dependent, which enabled us to address (i) the plateau issue on the evaluation of the FC and (ii) the non-Markovian behavior of the liquid-solid friction. For (i) we found that the Green-Kubo integral of the random force autocorrelation has actually a plateau for the finite-sized systems and the plateau value has a clear hydrodynamic meaning. For (ii) we evaluated the frequency-/time-dependent FC from equilibrium molecular dynamics simulation data for a Lennard-Jones (LJ) liquid and for water under different wall confinements and temperatures, without ambiguity due to the simulation system size. We showed that the Maxwell viscoelastic model is a fair approximation for the FC of LJ liquid on a fcc wall, and similarly a model with a single relaxation time can be applied for water on fcc and graphene walls at a high temperature, but a model with more than one time scale is required to describe the FC of supercooled water on both walls. Our theoretical framework opens the way to explore the frequency dependent FC for a wide range of complex liquids by non-demanding atomistic simulations, whose system size may be small. This work was financially supported by JSPS KAK-ENHI Grant Nos. 18K03929 and 18K03978, and by the ANR, Project ANR-16-CE06-0004-01 NECtAR. YY was also supported by JST CREST Grant No. JPMJCR18I1, Japan. LJ was supported by the Institut Universitaire de France.