Sound waves, diffusive transport, and wall slip in nanoconfined compressible fluids

Although continuum theories have been proven quite robust to describe confined fluid flow at molecular length scales, molecular dynamics (MD) simulations reveal mechanistic insights into the interfacial dissipation processes. Most MD simulations of confined fluids have used setups in which the lateral box size is not much larger than the gap height, thus breaking thin-film assumptions usually employed in continuum simulations. Here, we explicitly probe the long wavelength hydrodynamic correlations in confined simple fluids with MD and compare to gap-averaged continuum theories as typically applied in e.g. lubrication. Relaxation times obtained from equilibrium fluctuations interpolate between the theoretical limits from bulk hydrodynamics and continuum formulations with increasing wavelength. We show how to exploit this characteristic transition to measure viscosity and slip length in confined systems simultaneously from equilibrium MD simulations. Moreover, the gap-averaged theory describes a geometry-induced dispersion relation that leads to overdamped sound relaxation at large wavelengths, which is confirmed by our MD simulations. Our results add to the understanding of transport processes under strong confinement and might be of technological relevance for the design of nanofluidic devices due to the recent progress in fabrication methods.


I. INTRODUCTION
The field of nanofluidics has emerged from microfluidics within the last 20 years due to progress in nanofabrication, characterization and simulation tools [1][2][3].The reduction of characteristic length scales led to new phenomena that occur with increasing surface-to-volume ratios.A primary application of nanofluidic systems is single molecule selectivity in biotechnological lab-on-a-chip devices [4,5].The ultimate lower length scale limit of continuum theories has been determined to be 1 nm, and indeed, many effects occur on length scales well above this limit [3].However, recent advances in the fabrication of devices with sub-nanometer confinement [6][7][8] require theories that consider the discrete nature of molecules, and particularly, interactions between fluid molecules and the electronic structure of the confining walls become relevant [9,10].
The transition from continuum to molecular descriptions of confined fluids is closely related to fluid structuring.Molecular arrangement into distinct layers occurs due to geometric confinement and is strongly influenced by the interaction between the fluid and wall atoms [11][12][13].Molecular dynamics (MD) simulations have been widely employed to study the effect of fluid structure on transport properties [14][15][16][17][18] and the fluid-wall boundary condition [19][20][21].While MD simulations are a valuable tool to study confined system, Travis et al. [22] showed that classical continuum theories remain valid down to about five molecular diameters in simple fluids.Therefore, dimensionality-reduced continuum formulations are often employed on these scales, and combined with MD parameterizations, e. g. for wall slip [23].
In many continuum simulations of dense confined fluids, compressibility effects are neglected, which appears reasonable at first glance.Sound speeds in liquids are on the order of 10 3 m/s, and thus plays a minor role in momentum transport at hydrodynamic length and time scales [24].However, in the presence of walls, friction strongly influences the dynamics of sound waves, as was first pointed out by Ramaswamy and Mazenko [25] in a phenomenological approach for adsorbed layers on substrates.In their work, an unspecified friction term was used to interpolate between compressible hydrodynamics and Fickian diffusion.The effect of overdamped sound was later related to the (negative) algebraic long-time tail of the velocity autocorrelation function of a suspended particle in a confined system observed in lattice Boltzmann simulations [26,27].Although the diffusive sound modes govern the velocity of the suspended particle at long times, they do not contribute to the diffusion coefficient.A rigorous mathematical treatment of the problem was subsequently given in a series of works by Felderhof for a single plane wall [28], two parallel plane walls [29], and circular geometries [30].
It is common practice to extract transport coefficients from the correlations of equilibrium fluctuations, e. g. through the Green-Kubo approach [31,32], or by means of a direct fit to the autocorrelation functions of hydro-dynamic variables [33,34].
Here, we briefly recap bulk hydrodynamic theory in order to motivate a similar approach for confined fluids, see Refs.[35,36] for more details.We start from the Navier-Stokes equations that describe mass, momentum, and energy balance in terms of the density ρ( r), momentum j( r) and energy e( r) fields.Here and in the following, r = (x, y, z) denotes a spatial coordinate.Given the trajectory { r i (t), v i (t)} of N particles, we compute the mass density field as and the corresponding Fourier transform is given by Similarly, we obtain for the Fourier coefficients of the momentum density For small deviations from equilibrium δρ( k, t) = ρ( k, t) − ρ( k, t) , the linearized hydrodynamic equations can be solved.For the sake of brevity, we omit the δ-notation from now on.The time evolution is typically solved through a Laplace transform, and one finds that the longitudinal and transverse momentum modes decouple.The normalized time autocorrelation function of the Fourier coefficients of transverse momentum decays exponentially C ⊥ (k, t) ≡ j * ⊥ (k, 0)j ⊥ (k, t) j * ⊥ (k, 0)j ⊥ (k, 0) = exp(−νk 2 t), (4) where j ⊥ = j − j ˆ k and j = j • ˆ k with ˆ k = k/| k| are the momentum fluxes perpendicular and parallel to the wavevector k, respectively.The angular brackets denote an average over initial conditions, and the star is the complex conjugate.The decay rate is proportional to the kinematic viscosity ν = η/ρ, where η is the conventional shear viscosity.The kinematic viscosity has units of a diffusion constant, and can be regarded as such for the diffusion of momentum.
For the inverse Laplace transform of the longitudinal modes, one usually uses a second order approximation in k = | k| (instead of finding the exact roots of a cubic equation) [37].The normalized longitudinal momentum autocorrelation function has propagating sound modes and the decay rate is determined by viscous and conduc-tive effects with sound attenuation coefficient Γ = (γ − 1)D T /2 + ν L /2 and adiabatic speed of sound c s .Here, γ = c P /c V denotes the ratio of volume-specific isobaric and isochoric heat capacities, D T = κ T /c P is the thermal diffusivity with thermal conductivity κ T , and ν L = (4η/3 + ζ)/ρ is the kinematic longitudinal viscosity with bulk viscosity ζ.Finally, the normalized density autocorrelation function is given by where the first and second term describe the Rayleigh and Brillouin process, respectively.Equation.( 6) has probably received most attention, since its power spectrumthe dynamic structure factor S(k, ω)-is experimentally accessible, for instance through light scattering [38].The set of correlation functions fully describes the dynamics of statistically independent fluctuating quantities.The fluctuations of density are assumed to occur on time scales that do not allow the exchange of heat.Hence, the adiabatic compressibility governs the velocity of sound (D T k c s ) [35].In this paper, we investigate the time correlations of confined systems.Gutkowicz-Krusin and Procaccia [39,40] have derived the dynamic structure factor for confined fluids, considering both momentum and energy transport at the fluid-wall interface.However, their approach considers only no-slip or perfect slip boundary conditions.Bocquet and Barrat [41,42] derived momentum time correlation functions for confined systems, and related these to the hydrodynamic boundary conditions with partial-slip.Yet, they considered only transverse fluxes perpendicular to the walls, averaged over the lateral dimensions of the system.The effect of the lateral periodic box size in MD simulations of confined systems, i. e. the admissible wavelengths of longitudinal and inplane transverse modes, have received little attention so far [43].
Here, we provide a derivation of time correlation functions for height-averaged balance equations-similar to the continuum methods commonly applied in lubrication [44]-under isothermal conditions.Thus, we consider only momentum transport at the interface with partialslip boundary conditions, which allows us to extract effective transport coefficients as well as the slip length from equilibrium MD simulations.Under nanometric confinement, we are able to probe the transition to overdamped sound modes at long wavelengths.

A. Dimensionality reduction
Continuum formulations of thin-film flows typically employ an average over the gap height to reduce the dimensionality of the problem.We have shown in a previous work [45], how this can be formally achieved without making a priori assumptions about the constitutive behavior, which we briefly recap here.We average the hydrodynamic balance laws over the gap height where q ≡ q( r, t) denotes the vector of conserved densities, and f i ≡ f i ( r, t) are the corresponding fluxes in Cartesian direction i.Note that here and in the following, bold symbols (e.g.q, f ) indicate vectors of arbitrary length representing a collection of state variables (or derived quantities) while arrows (e.q.r, k) indicate Cartesian 3-vectors.Formally, this integration can be performed for channels with surfaces moving lateral to each other and having surface topography, such as in lubrication.Hence, the integration limits h 1 (x, y, t) and h 2 (x, y, t) depend on the lateral Cartesian coordinates and on time.We denote the gap height with For the integral on the l.h.s. of Eq. ( 7) and the first two terms on the r.h.s.Leibniz rule for differentiation under the integral sign applies, and after a few steps [45] one arrives at a dimensionality-reduced form of the balance equations where overbars denote height-averages φ = 1 h h2 h1 φ dz, and s acts as a source term.Due to the structure of Eq. (8), where dominant diffusive fluxes (e. g. shear stresses) are lumped into the source term, explicit numerical schemes for hyperbolic balance equations with source terms, have been proven successful to solve lubrication problems [45].
In the most general case, the source term contains flux boundary conditions in the averaging direction, as well as terms which depend on the topography and movement of the upper and lower wall.Here, considering only flat channels without shearing, the source term simplifies to The source term vanishes when there is no flux across the boundary, i. e. for impenetrable, perfectly insulating, and slippery walls, where Eq. ( 8) describes a two-dimensional fluid.For nonzero source terms, additional dissipation relative to a laboratory system is added.In the following, we focus on systems with impermeable walls under isothermal conditions.In order to solve Eq. ( 8), we need explicit expressions of the relevant fluxes as a function of the conserved variables, f i (q) (or their averaged versions), i. e. constitutive relations.

B. Hydrodynamic correlations in confined systems
In the following, we focus on isothermal conditions where the density vector q = (ρ( r, t), j( r, t)) does not contain the energy density.correlation functions of density and longitudinal momentum.For bulk fluids, this special case is recovered by setting γ = 1, which renders adiabatic (c s ) and isothermal sound speed (c T ) equal.Furthermore, we neglect nonlinear convective terms, which is justified by the thin-film assumption [46].The flux in Cartesian direction i ∈ [x, y] and the source term are then given by fi respectively, where δ ij is the Kronecker symbol.Here, τ ij denotes the components of the viscous stress tensor, which for a Newtonian fluid in three dimensions reads where η and ζ are the coefficients of shear and bulk viscosity respectively, u is the velocity field, and 1 is the 3×3 unit matrix.The pressure p is given by a barotropic equation of state (EOS) p(ρ) = c 2 T ρ, which was chosen to retain the linearity of the problem.
In nanoscale geometries, deviations from the no-slip boundary conditions become relevant [47][48][49].Therefore, we consider Navier slip boundary conditions [50] with a uniform slip length b both at the top and bottom surface.The slip length is the virtual distance from the fluidwall interface at which the fluid velocity reaches the wall velocity, if linearly extrapolated.With the assumption that density does not vary across the gap, we obtain the height-averaged fluxes and the source term Here, κ renormalizes the actual gap height h for a system with slip to an effective gap height h eff = h √ κ for an equivalent system without slip.The expression for κ can be obtained by shifting a parabolic Poiseuille velocity profile u(z) to no slip boundary conditions while maintaining the same average flux as in the slip-profile, i. e.
where u * (z) denotes the shifted profile.A brief derivation of κ for different slip lengths at the top and bottom wall is given in Appendix A. In the scope of this work, we only deal with the case of identical slip length b at both walls, for which we obtain and which is illustrated in Fig. 1.
Wall slip leads to an effective gap height h eff = h √ κ.The parameter κ is found by equating the z-average of the original velocity profile u(z) with that of u * (z), which is u shifted to yield zero slip.
To arrive at a general solution to Eq. ( 8), we express the densities of conserved variables as a series of normal modes q( r, t) = q( k, t)e i k• r , where k and r are twodimensional vectors in the plane of the confined region.Note, that for convenience, we keep the notation introduced above, but all vector dimensions are reduced by one due to the average.Since all field variables are averaged over the gap height, we drop overbars from now on for the sake of brevity.We obtain an ordinary differential equation for the Fourier coefficients with the hydrodynamic matrix The hydrodynamic matrix H is diagonalized with eigenvalues corresponding to transverse modes, and corresponding to longitudinal modes.Here, the isothermal speed of sound s T (k) = c 2 T − (τ k) −2 follows a confinement-induced dispersion relation with τ = (Re µ ) −1 which for small wavenumbers (hk 1) reads Surprisingly, the dispersion takes effect at large wavelengths λ = 2π/k, where the second term in the discriminant becomes important, and the speed of sound deviates from its bulk counterpart given by the EOS.As long as s T is a real number, longitudinal modes show underdamped oscillations.However, there is a transition from underdamped to overdamped behavior at a critical wavelength where s T becomes imaginary and the eigenvalues for longitudinal modes µ become real.In Fig. 2, we illustrate the wavelength dependence of the effective speed of sound.Note, that we plot absolute values for s T normalized by kτ as a function of the wavenumber normalized by k crit .Similar expressions can be derived for axisymmetric flow through circular channels with radius R, where the only two eigenvalues are given by with κ = 1 + 4b/R.The definition of derived quantities such as λ crit change accordingly.
Consequently, modes with wavelengths larger than the critical wavelength λ crit cannot propagate through the channel and are dissipated within a finite relaxation time.This is rather counterintuitive, since in bulk hydrodynamics wavelength-independent transport coefficients are assumed at sufficiently long wavelengths.
Assuming arbitrary initial conditions q(k, 0) = (ρ(k, 0), j (k, 0), j⊥ (k, 0)) (with the limitation that perturbations out of equilibrium are small in order to preserve linearity), we solve Eq. ( 16) for the time evolution Dispersion relation for a height-averaged continuum formulation of confined fluids.The solid line describes the angular frequency using the magnitude of the phase velocity, given by Eq. ( 20).The dashed line describes the bulk reference.The wavenumber is normalized by kcrit = 6ν/h 2 κcT and the angular frequency is normalized by the characteristic relaxation time of longitudinal modes in the underdamped limit lim k→k crit τ = h 2 κ/6ν.The critical transition from underdamped to overdamped dynamics occurs with diverging group velocity.
of the real part of the conserved variables, with characteristic relaxation times Note, that these expressions remain valid when s T becomes purely imaginary and the trigonometric functions turn into their hyperbolic counterparts.We give the solution of the imaginary parts explicitly in Appendix B for completeness.Figures 3a and b show Eq. ( 23a) and (23b), respectively, as well as their dependence on wavelength in the case of underdamped oscillations (λ < λ crit ).Since λ crit ∝ c T h 2 /ν, and c T /ν ∼ O(1) for most dense fluids, the critical wavelength is on the order of the magnitude of the squared gap height.In the overdamped case (λ > λ crit ), s T is an imaginary number, and therefore, the behavior of density and longitudinal momentum modes changes fundamentally as shown in Fig. 3c-d.We observe that density perturbations do not decay, i. e. their relaxation time diverges.Furthermore, one can show that the decay rate scales with k 2 , similar to bulk fluids (see Appendix C).However, in contrast to the bulk, sound relaxation times converge to a finite value, identical to that of transverse modes.We want to highlight, that the functional form of Eq. ( 23) is equivalent to the autocorrelation function of equilibrium fluctuations in the bulk, but with different characteristic time scales, as described above.In particular, relaxation times of longitudinal and transverse modes become wavelength-independent in the long wavelength limit, and a diverging group velocity leads to a transition to overdamped sound relaxation.We explicitly probed this transition using equilibrium molecular dynamics simulations of confined fluids.Therefore, we write generic correlation functions of the form where β = 1/ωτ , and use frequencies and decay rates as fit parameters to the autocorrelation functions obtained from MD.Note that for bulk systems, since β ∼ O(k), the sine-term is often neglected, but some authors [51][52][53] explicitly include it.As highlighted in the beginning of this section, our derivation is for non-fluctuating isothermal conditions with γ = 1.However, the numerical tests with MD naturally contain thermal fluctuations at macroscopic constant temperature, which is why we included the thermal part in Eq. ( 25c).An overview of the theoretical expressions for ω and τ for different systems and conditions is given in Tab.I.

III. MOLECULAR DYNAMICS SIMULATIONS
In the previous section, we derived expressions for the Fourier coefficients of conserved variables in confined fluids in the hydrodynamic limit.To scrutinize our predictions, we performed a brute force test of our findings using molecular dynamics (MD) simulations of a simple fluid, confined in a nanometer-sized channel.All MD simulations were carried out with LAMMPS [54], and we have used a supercritical Lennard-Jones fluid at the state point T = 2.0 /k B and ρ = 0.452 σ −3 for all results shown in this paper.Yet, our findings are not limited to the supercritical state, which we showed in a related work [55].Hence, the interaction potential [56] between fluid atoms is given by where r ij = | r i − r j | is the distance between particle i and j.The interatomic potential is shifted to zero at a cutoff radius of r c = 2.5σ.All simulations were performed with a time step ∆t = 0.0025 mσ 2 / .We performed bulk reference simulations in fully periodic boxes, as well as simulations in confined system, where the periodicity is broken in the direction normal to the walls, which are modeled as explicit rigid atoms.To probe the long wavelength limit, we used simulation boxes with high aspect ratios, i. e. where one box length in the direction parallel to the walls is much larger than both the gap height and the remaining in-plane dimension.We sampled the trajectories in the microcanonical ensemble for 2 × 10 7 time steps after an initial equi-libration of the system in the canonical ensemble, and recorded every 2000th step for post-processing.

A. Confined setup
For the confined fluid simulations, we model the walls as two rigid atomic layers arranged in an fcc lattice with lattice constant a = 1.2σ f = 1.2σ, which corresponds to a wall density of ρ w = 2.31σ −3 .The wall density significantly affects fluid-wall commensurability, hence slip [19], and the chosen parameter lies within the range of similar studies [19,57,58].The [1 1 2] and [ 1 1 0] directions are taken in x-and y-direction, respectively, such that the closest-packed {111}-surface is in contact with the fluid.In x-and y-direction, we keep periodic boundary conditions, hence the lateral box sizes are always chosen to be multiples of 3a/ √ 2 and a/ √ 2 respectively.Fluidwall interactions are governed by a Lennard-Jones potential.We use the Lorentz mixing rule for the length scales σ wf = (σ w + σ f )/2, with σ w = 0.75σ.We model different wetting behavior at the interface by scaling the fluid-wall interaction with respect to the fluid's interaction energy, i. e. wf = α f , where the parameter α can be related to the contact angle [59].
Input parameters for the confined system simulations are gap height, lateral box sizes, and number density, which defines the number of atoms N in a homogeneous system.Note, that the inner wall layers are placed at a distance h + σ apart to account for a thin depletion zone at the interface.In the case of pronounced layering effects, fixing the gap height can lead to deviations from the target density.Yet, we stick to this pragmatic approach, which is particularly useful when comparing atomistic results with continuum predictions.For the gap heights and wetting properties covered within this study, we did not observe strong deviations from the target density in the center of the fluid film.The confined MD setup is illustrated in Fig. 4 LJ fluid

B. Autocorrelation functions from equilibrium fluctuations
For both bulk and confined systems, we compute the autocorrelation functions for the Fourier coefficients of mass and momentum density.We compute the average in Eq. ( 4)-( 6) from a single time series of the dynamical variables, which is equivalent to the convolution of that time series with itself.We leverage the computational efficiency of the Fast Fourier Transform (FFT) and use the Wiener-Khinchin theorem to compute the integral [60].We chose k = (k n , 0, 0) in the direction of the longest dimension of our simulation box, where k n is a discrete wavevector corresponding to waves that fit into this periodic dimension, i. e. k n = 2πn/L x n = 1, 2, 3, . ... Hence, momentum autocorrelation functions are calculated for the transverse ( j ⊥ k) and longitudinal ( j k) direction.
To rule out possible size or shape effects in this setup, we varied the lateral box sizes independently and compared transverse momentum autocorrelation functions computed from either of the lateral velocity components, which led to indistinguishable results.Hence, computing transport coefficients from correlation functions of collective variables in boxes with high aspect ratio seems to be no issue, in contrast to e. g. the calculation of selfdiffusion coefficients [61,62].

C. Non-equilibrium simulations of slip
We performed non-equilibrium MD simulations with a similar setup as described in Sec.III A. Instead of sampling the equilibrium state, we sheared the walls at a constant shear rate, to obtain reference values for the slip length.Therefore, we used smaller box sizes in the shearing direction (L x = 29.40σ)resulting in 2818 fluid atoms.We sheared the upper and lower wall at constant velocity ±u 0 /2, respectively, and sampled the Couette profile after an initial startup period until a total sliding distance of 100L x is reached.
Shearing at a fixed gap height with rigid wall atoms may lead to shear localization [63] or diverging slip lengths [64].The choice of our simplified setup is mainly to ensure better comparison with equilibrium simulations, where these effects do not play a role.Here, we did not observe such phenomena at the densities and shear rates considered.Furthermore, the maximum applied shear rate u 0 /h in our reference simulations is lower than 0.01 /mσ 2 , where the shear rate dependence is expected to be low, even for strong fluid wall interaction [65].

D. Green-Kubo simulations
We performed bulk equilibrium simulations in cubic boxes with 1000 atoms to obtain reference viscosities using the Green-Kubo approach.The GK integrals for the shear and bulk viscosity read respectively, where τ ij are the components of the deviatoric stress tensor and δp(t) = tr [σ(t)] /3− p with stress tensor σ(t) and the equilibrium average of the hydrostatic pressure p .The value of the integrals in Eq. ( 27a) and (27b) converges after approximately 500∆t, for correlation functions computed from a trajectory with 4 × 10 5 time steps.The viscosities have been computed from the integrals of 20 equivalent replica simulations.
IV. RESULTS

A. Autocorrelation functions
We first computed the autocorrelation functions from a bulk reference simulation without walls.The simulation box had dimensions 941.2 × 14.7 × 14.7 σ 3 , which corresponds to 92 001 atoms.Confined system simulations were performed in a similar box with 90 189 fluid atoms and 43 520 solid atoms in each of the two walls.In this section, we show simulation results with a wall-fluid interaction parameter α = 0.75.Note, that typical values for σ are in the range of a few Ångströms, and therefore, our systems have gap heights of a few nanometers.We obtained characteristic times for the decay and propagation of equilibrium fluctuations by a least-square fit to the generic expressions, Eq. ( 25), for discrete wavenumbers.We used MD data in a time interval of 250 mσ 2 / for the fit, and for the confined system all wavelengths considered at this point are in the underdamped regime.In all cases, we used rates (e.g. 1/τ ⊥ ) and frequencies(ω) as fitting parameters, but plot the inverse, i. e. the corresponding characteristic times.
Transverse momentum correlations are independent of the longitudinal modes, and the kinematic shear viscosity ν is the corresponding transport coefficient.Figures 5ad show the autocorrelation functions of the real part of j y (k, t) for four wavelengths (blue solid lines) and the corresponding fit (black dashed lines), where the shear relaxation time τ ⊥ increases with wavelength.
From the transverse momentum autocorrelation functions of the confined fluid in Fig. 5a-d it becomes evident, that the wavelength dependence in the decay rate is lost in the range of the presented wavelengths.This is qualitatively in agreement with the prediction in Eq. (23).
The relaxation times of transverse momentum fluctuations are much shorter than in the bulk.
The autocorrelation functions of longitudinal momentum fluctuations obtained from MD simulations as well as the corresponding fits are shown in Fig. 5e-h.Different relaxation behavior between the bulk and the confined systems can be seen in the longitudinal direction as well, with stronger damping in the system with walls.For the shorter wavelengths, the positions of the local minima and maxima in the autocorrelation functions coincide with those of the bulk fluid.For the largest wavelength in Fig. 5e, we observe that the first minimum is shifted to shorter times compared to the bulk, indicating a frequency shift with increasing wavelengths.
The density autocorrelation function Eq. ( 6) contains the previously obtained longitudinal momentum autocorrelation function weighted by 1/γ.Therefore, to reduce the number of fitting parameters, we fix the sound attenuation rate and sound frequency, which leaves only the decay rate of thermally induced density fluctuations 1/D T k 2 and the heat capacity ratio γ to be determined.The results are shown in Fig. 5i-l, where we make similar observations of faster relaxation with weaker wavelengthdependence as in the bulk fluid.

B. Effective time scales under confinement
The correlations of equilibrium fluctuations in confined fluids clearly deviate from their bulk counterparts, as shown in the previous section, but Fig. 5 highlighted only four selected wavelengths.In the following, we systematically investigate the transition of characteristic time scales as wavelength increases.In particular, we test whether the isothermal, height-averaged theory described in Sec.II A adequately describes the relaxation times and frequencies obtained from MD.Therefore, we computed all quantities appearing in Eq. (24b), (24b), and (20) in separate non-equilibrium (see Sec. III C) and equilibrium MD simulations (see Sec. III D), and compare the theoretical predictions with the corresponding fit parameters from the autocorrelation functions.Since, thermal effects are not considered here, but are naturally included in the MD simulations, we obtained additional parameters (C P , C V , κ T ) from the NIST thermophysical database for Argon [66].A summary of material parameters at the supercritical state point is given in Tab.II.
Figure .6a shows the shear attenuation time τ ⊥ over the wavelength as obtained from both bulk and confined MD simulations.The symbols represent an average over the two transverse directions (j y (k, t) and j z (k, t), real and imaginary parts respectively) in the bulk, and over the in-plane transverse direction (from Re[j y (k, t)] and Im[j y (k, t)]) in the confined system.We do not investigate correlations of the j z -Fourier coefficients, since we assumed laminar flow to motivate the height-averaged balance equations in Sec.II A.
Shear relaxation times in the bulk scale with the square ‡ from NIST [66] of the wavelength as a consequence of momentum conservation.The dash-dotted line in Fig. 6a illustrates the prediction for the bulk fluid with kinematic viscosity as the constant of proportionality.As expected, the shear relaxation time in confined systems deviates from the bulk with increasing wavelength, and converges to a constant value of approximately 30 mσ 2 / for wavelengths larger than 120σ.The prediction for the confined fluid based on Eq. ( 24a) is shown as a dashed line, and is in excellent agreement with the MD data.Similar behavior can be observed for the sound attenuation time τ in Fig. 6b.Here, symbols for both bulk and confined configuration correspond to an average over the fitted relaxation times from the real and imaginary part of the longitudinal momentum correlations.The bulk relaxation times clearly follow a quadratic scaling relation, and the dash-dotted line illustrates the prediction with shear-attenuation coefficient Γ given in Tab.II.For the confined fluid, the transition to wavelength-independent relaxation times is similar to the in-plane shear relaxation.Here, the characteristic time for the decay of sound modes converges to a value approximately twice as large as the shear relaxation time in the long-wavelength limit, and to the same quadratic scaling as in the bulk in the limit of short wavelengths.The prediction based on the isothermal theory, Eq. (24b), is shown as a dashed line.In the long wavelength limit, relaxation times from MD and the predictions agree, but the isothermal theory underestimates sound attenuation at shorter wavelengths.Therefore, we plot a modified prediction as a dotted line, which replaces the first term of Eq. (24b) with the nonisothermal sound attenuation from the bulk theory, which describes the MD data slightly better.
For completeness, we show the relaxation behavior associated with the Rayleigh process, i. e. thermal relaxation times obtained from the density autocorrelation functions, in Fig. 6c.The dash-dotted line describes The first column (a-d) corresponds to transverse momentum ( j ⊥ k), the second column (e-h) corresponds to longitudinal momentum ( j k), and the third column (i-l) corresponds to density correlations.Panels within the same row share the same wavelength.Solid lines show results from molecular dynamics (MD) calculations, and dashed lines are the corresponding fits.For simulations of the confined system, all wavelengths shown here are in the underdamped regime (470.6σ< λcrit).
again the theoretical expectation based on the bulk theory, where the thermal diffusivity D T is entirely based on data from NIST.Thermal relaxation times from bulk and confined MD simulations both scale with the square of the wavelength and are indistinguishable from one another in the case of rigid walls.We show additional data obtained from simulations in shorter boxes with vibrating wall atoms with open symbols, which indicates a similar transition to wavelength-independent relaxation times if energy transport through the walls is allowed.
Oscillatory behavior of density and longitudinal momentum correlation encodes sound propagating properties of the considered fluid systems.In Fig. 7a, we plot the sound period T = 2π/ω over the wavelength for both bulk and confined fluids, obtained from fits to the theoretical expression for the momentum autocorrelation function.For illustrative purposes, we chose stronger wall-fluid interactions with α = 1.5, as compared to α = 0.75 in Fig. 5 and 6, since deviation from bulk behavior becomes more pronounced.As expected, the sound period obtained from the bulk MD system scales linearly with the wavelength and matches the prediction (dash-dotted line) based on the adiabatic sound speed from Tab. II.
For small wavelengths, the sound period in the confined fluid follows approximately the same linear relation as the bulk, with a slight positive shift to longer periods-or lower sound speeds.With increasing wavelength, deviation from bulk behavior becomes more pronounced.This effect becomes clearer when we plot the phase velocity directly over the wavelength in Fig. 7b and compare to the dispersion relation Eq. ( 20) derived in Sec.II A. As expected, the phase velocity is constant for bulk fluids in the considered range of wavelengths.
Simulations of the confined system reveal constant speed of sound only for small wavelengths up to approximately 120σ.For larger wavelengths, the sound speed decreases, and we found good agreement with the dispersion relation Eq. ( 20).Using bulk reference data from Tab. II and Eq. ( 21), we predict a critical wavelength λ crit = 675σ.

C. Critical damping and overdamped regime
As we have described in Sec.II A, the functional form of density and momentum correlations in confined systems change fundamentally at the critical wavelength, when s T becomes imaginary and the trigonometric functions turn into their hyperbolic counterparts.In the overdamped regime, the decay of sound modes is governed by two distinct relaxation times for density and momentum correlations.Therefore, we seek approximations to Eq. (25b) and (25c) which are more suitable for fitting to the available MD data at large wavelengths.We show in Appendix C how to arrive at effective time correlation functions in the overdamped regime, and distinguish between τ j and τ ρ for the relaxation times of longitudinal momentum and density fluctuations, respectively.Thermal relaxation is slow at large wavelengths, such that we can assume exp(−D T k 2 t) ≈ 1 in the considered time interval.Hence, we fitted the time autocorrelation functions of density and longitudinal momentum in the overdamped regime to the following expressions The critical wavelength for our confined MD simulation with gap height h = 14.7σ and α = 1.5 is approximately 675σ-more than one half of the full box length (941.2 σ)-such that only the largest wavelength remains available for analysis of the overdamped relaxation.Therefore, we further increased the box length to L x = 1411.8σ,resulting in two more data points in the overdamped regime.Figure 8 illustrates the split into two distinct relaxation times for density and longitudinal momentum, which were formerly identical in the underdamped region.The dashed lines illustrate the theoretical prediction based on the eigenvalues corresponding to longitudinal and density modes Eq. (19).The data points available from our MD simulations converge towards the theoretical predictions, and the corresponding autocorrelation functions are shown in the inset.

D. Continuum simulations
To further illustrate the suppression of sound waves in the overdamped regime, we performed continuum simulations of confined fluids with a finite volume implementation of Eq. ( 8).The details of the implementation can be found in Ref. [45].We studied the dynamics of density perturbations through explicit time integration in one-dimensional slit channels with length L x = 1 µm and varying gap heights between 5 and 50 nm.We employed no-slip conditions at both channel walls and chose material properties for supercritical argon, similar to our MD simulations.The explicit time step was ∆t = 500 fs and 1024 grid cells discretized the domain.We initialized the system with an equilibrium density ρ 0 and added a Gaussian-shaped density perturbation δρ(x) centered in the middle of the channel with standard deviation σ = L x /40.The simulation time was short enough to avoid effects from the periodic boundary conditions.
Figure 9 illustrates the two limiting behaviors of sound propagation in the underdamped regime and diffusion in the overdamped regime.The gap height in Fig. 9a-b is ten times larger than in Fig. 9c-d, while all other parameters are kept constant.Note, that in both cases there is no single well-defined wavelength, but we take the width of the wave package as the dominant contribution and use it to distinguish between underdamped and overdamped behavior.
In the underdamped case, the initial Gaussian density profile separates into two wavelets propagating to Transition from underdamped to overdamped dynamics.The last four sound relaxation times in the underdamped regime are shown as green diamonds.After reaching the critical wavelength, we obtain two real eigenvalues from Eq. ( 19).One leads to a converging relaxation time for momentum perturbations (lower branch of the dashed line), the other one to a diverging relaxation time for density perturbations (upper branch of the dashed line).For α = 1.5 and box lengths of 941.2σ and 1411.8σ,we characterize three modes in the overdamped regime with our MD calculations.The fits to the simplified (exponential) autocorrelation functions are shown in the inset, and the resulting relaxation times of density and momentum modes are shown as upward and downward pointing triangles, respectively.
the left and right.Here, we show only the right half of the symmetric profiles.Figure 9b shows the propagated distance over time.Hence, the slope tells us the speed of sound, which unsurprisingly yields c T , the isothermal sound speed as defined by the initial choice of an equation of state.
In the overdamped case, the same initial density perturbation decay is purely diffusive.Here, we measured the variance σ 2 of the density distribution as shown in Fig. 9d.The linear relationship clearly illustrates the diffusive transport, and the diffusion coefficient is given by D ρ = σ σ(t).The slope is unity in the normalized representation, which leads to

E. Effective gap height and slip length
Transverse momentum fluctuations do not propagate, and are therefore unaffected by the transition at λ crit .Due to their simple functional form and the excellent agreement of predicted and measured relaxation times 9. Time evolution of a small, Gaussian-shaped density perturbation in a slit geometry.The gap height in panel a is large enough to accommodate sound waves.Panel b shows the propagated distance ∆x over time t, rescaled by the critical wavenumber kcrit and the long-wavelength sound attenuation time τ = lim k→k crit τ , respectively.Therefore, we recover the sound speed cT = 1/τ kcrit from the slope of the curve.In panel c, the gap height is one order of magnitude smaller, which leads to a purely diffusive behavior.Panel d shows the variance σ 2 over time, rescaled by the squared critical wavenumber and sound attenuation time, respectively.Hence, the slope is proportional to the corresponding diffusion coefficient Dρ = h 2 c 2 T /12ν.Only the r.h.s. of the symmetric profiles is shown in a and b.
in Fig. 6a, we propose the following method to quantify viscosity and slip in confined systems using long wavelength correlations.Instead of predicting the transition from quadratic scaling to constant relaxation time using bulk and interfacial properties as in Sec.IV B, we fit Eq. (24a) to data obtained from MD.Thus, the characteristic transition enables us to obtain both viscosity and slip length as fitting parameters.The results for such a fitting procedure for the slip length is shown in Fig. 10a for various fluid-wall interaction energies.
The highest shear relaxation time can be observed when a purely repulsive Weeks-Chandler-Anderson (WCA) potential [67] is used for the wall-fluid interaction.The relaxation time decreases with increasing strength of attractive forces between fluid and wall atoms.
With increasing wall fluid interaction, from purely repulsive to α = 1.5, the effective viscosity increased by about 8.1 %.Compared to the bulk, the kinematic shear viscosity increased by 13.0 %.From the effective gap height, we obtain the slip length b from Eq. ( 15), and the results are shown in Fig. 10b.For comparison, we performed non-equilibrium MD simulations of identical (but shorter) systems using moving walls.We extracted the velocity profiles from these shear simulations by computing time averages over slices along the gap coordinate, and calculated the slip length from fits to the linear Couette profiles.We obtained similar slip lengths for both methods, however, with increasing wall-fluid interaction, the slip lengths from non-equilibrium runs tend to zero slip, whereas the equilibrium slip lengths start to saturate at α ≥ 1.0.

V. DISCUSSION
Onsager's principle [68], the fact that small fluctuations around an equilibrium state diffuse and propagate in the same way as large non-equilibrium perturbations, is the foundation of several methods to obtain transport coefficients from time correlations of collective variables at equilibrium, with Green-Kubo methods [31,32] be-ing the most prominent ones.In this work, we exploit this theoretical framework to study hydrodynamic correlations in confined systems.In contrast to the seminal work by Bocquet and Barrat [41,42], we considered averages over the confining dimension and explicitly probed long wavelengths using elongated simulation boxes.
The autocorrelation functions of density and momentum fluctuations in confined fluids have the same functional form as in the bulk, but transport coefficients strongly deviate from the bulk expressions.In a previous work, Porcheron and Schoen [53] came to the same conclusion, namely that the functional form of the correlation function does not depend on confinement, but the transport coefficients therein do.However, they considered only pure slip boundary conditions, such that momentum fluxes at the walls disappear.In this case, averaging the hydrodynamic equations over the gap height leads to a "true" two-dimensional fluid with vanishing source term in Eq. (8).Hence, the transition to finite relaxation times and the overdamped sound regime do not emerge, and are difficult to observe in MD simulations if small boxes are used, as shown in Figs. 6 and 7.Then, the interface only implicitly affects the relaxation constants, e. g. by viscosity changes due to commensurability effects between the gap height and the periodicity of fluid layering [53].A clear distinction between interfacial and fluid friction as recently suggested by Zhou et al. [69] is therefore difficult, when ordering effects reach far into the fluid region and confining dimensions are small.
Ordering effects are small in the presented MD simulations of confined LJ fluids, and therefore prediction of the spectral relaxation behavior based on bulk fluid properties works well in most situations.The prediction works particularly well for transverse modes, as shown in Fig. 6a, where the shear viscosity and slip length has been obtained from reference simulations in equivalent MD setups.In fact, fits to transverse momentum autocorrelation functions as proposed by Palmer [33] is an established method to calculate the viscosity at equilibrium [70] and has been recently extended by Cheng and Frenkel [34] to compute heat conductivity from density correlations.
Quantification of slip in MD simulations is of great importance for the understanding of confined fluid systems, e. g. for multiscale simulations of fluid transport [71] or lubrication [23].The characteristic transition from bulklike to wavelength-independent relaxation offers a new path to determine the slip length at equilibrium, as illustrated in Fig. 10.The proposed method is in the spirit of the early work by Palmer [33] and has the advantage that both interfacial slip and the shear viscosity of the fluid can be determined simultaneously.This can be seen as extension to a similar approach by Sokhan and Quirke [72], where viscosity is an input parameter, that has to be taken from literature, or is obtained from separate MD simulations.
Figure 10a highlights that even for purely repulsive fluid-wall interactions relaxation times eventually con-verge to a finite limit with increasing wavelength.Hence, the limit of infinite slip, where Galilean invariance is fully restored, is impossible to reach with MD simulations of interacting fluid and wall atoms.However, in the case of ultra-low interfacial friction, it requires large wavelengths to observe deviations from a bulk system which are usually not probed in conventional MD systems.
The prediction of the slip length based on the equilibrium simulations have been compared with results from non-equilibrium calculations in Fig. 10b.Slip lengths from both calculations decay exponentially with increasing wall fluid interaction.Equilibrium results seem to converge towards a finite slip length of approximately one atomic diameter, which is probably the closest one can get to a no-slip boundary condition in an atomistic system.This agrees with other equilibrium slip length measurements [73], suggesting that friction converges and leads to a small but finite slip length, even for strong fluid-wall interaction.Non-equilibrium measurements rely on linear extrapolation of the velocity profiles to obtain the slip length, which is very sensitive to the exact location of the fluid-wall interface.This might explain the deviation of the non-equilibrium results at higher interaction energies.For practical use of the proposed method, shorter box sizes on the order of a few gap heights might be sufficient, which drastically reduces the computational effort compared to the MD setup shown here.
Longitudinal momentum relaxation times show a similar behavior as in the transverse direction, as shown in Fig. 6b.Yet, in contrast to the transverse case, sound absorption depends on viscous and thermal effects.Therefore, predicted relaxation times for the bulk system based on literature data and our own simulation are slightly above the ones obtained from the autocorrelation functions.The isothermal theory for confined systems accurately describes the long wavelength relaxation, since the time required for heat to diffuse along the lateral direction is long enough to assume isothermal conditions.However, at short wavelengths, the theory overestimates relaxation times, but the prediction can be improved by considering thermal effects in the bulk contribution to the overall relaxation.Remaining deviations between our prediction and the data, might be due to viscosity enhancement in the confined system compared to the bulk due to ordering effects introduced by the walls, as discussed above.
The role of thermal effects in the presented MD results is highlighted in Fig. 6c.Relaxation times for the Rayleigh process contribution in the density autocorrelation function obtained from bulk and confined systems are indistinguishable, since rigid walls are used.In this context, rigid walls imply that there is no thermal coupling between fluid and wall atoms, since wall atoms do not vibrate, i. e. there is no heat flux across the interface.In analogy to wall slip described above, this represents an idealized situation with a perfect "thermal slip" condition, i. e. an interface with infinite Kapitza length [74].
Therefore, similar transition to finite relaxation times is expected for systems with vibrating walls, which is shown in Fig. 6c, but is not the main focus of this paper.
Instead, we focus on the long-wavelength behavior of propagating modes.The dispersion relation, Eq. ( 20), is dominated by the longitudinal relaxation time in the long wavelength limit and therefore independent of thermal effects.Hence, considering vibrating walls does not affect the existence of the transition to overdamped sound.As shown in Fig. 7, both sound period and accordingly the velocity of sound in confined fluids can be adequately described with Eq. (20).In bulk fluids, dispersion is usually observed at molecular lengths scales, where transport coefficients depend on wavelength, and short-term memory effects have to be considered within a generalized hydrodynamic theory [35].In confined systems, dispersion is an intrinsic property of the system, where both fluid and interfacial properties determine the lengthscale at which it becomes effective.
The predicted bifurcation into diverging density and converging momentum relaxation times in the overdamped regime has been confirmed by our simulations in very large systems, as shown in Fig. 8. Hence, long wavelength modes in confined systems are transported entirely by diffusion.The short-lived longitudinal momentum correlations can be interpreted as discrete "jumps" of a wave package, that bring about diffusive mass transport in hydrodynamic systems with additional dissipation due to walls, similar to Fickian diffusion for discrete particles.The size of typical MD simulations of confined fluids does not probe this limit, although it might be the dominant transport mechanism in real confined systems.For instance, Cheng and Giordano [48] measured slip lengths of 25 nm in a 50 nm thin channel for hexadecane on photoresist-coated glass.The critical wavelength for this system is on the order of micrometers, thus reaching into the frequency range of ultrasound applications.
The results of our continuum simulations shown in Fig. 9 illustrate the effect of overdamped sound in a non-equilibrium scenario.By reducing the gap height, we interpolate between sound propagation and diffusion, as initially suggested by Ramaswamy and Mazenko [25].The empirically determined diffusion coefficient agrees with the one reported in Ref. [75].Transition to overdamped behavior may also lead to arrest of an initially propagating wave package due to spread related to dissipation in the underdamped regime or due to spatially varying gap height or wall slip.The latter suggests further work on the role of roughness on the critical transition of sound modes in confined systems.

VI. CONCLUSION
In this work, we showed that correlations of the hydrodynamic conserved variables in confined fluids can be derived from an isothermal height-averaged description of continuum balance equations.The functional form of hydrodynamic correlation functions remains equivalent to the bulk, but characteristic time scales therein are affected by confinement.We focused on the lateral wavelength dependence of relaxation times and phase velocities.The continuum description predicts a transition to constant relaxation times of density and momentum perturbations, which we confirmed by upscaling MD simulations to the long wavelength limit.Furthermore, our theory contains a geometry-dependent dispersion relation for the speed of sound in the long wavelength limit, which leads to a transition from underdamped to overdamped dynamics, which is also evident from the MD simulations.Large MD boxes are required to probe the overdamped limit, but the transition can be on the order of the system size in highly confined fluids.Hence, diffusive contributions to lateral mass transport might be systematically neglected in finite systems.Finally, we proposed a new equilibrium method to calculate the slip length in confined MD systems based on our findings, which shows accurate results when compared to a non-equilibrium reference. and in Eq. (A1).Evaluating Eq. ( 14) eventually leads to Appendix B: Imaginary parts The imaginary parts of the solution to Eq. ( 16) are shown here for completeness.and we are now interested how this limit is reached.For the sake of brevity, we write a ≡ 6ν/h 2 κ and b ≡ −is T = (a/k) 2 − c 2 T with c T < a/k, i. e.
ρ(k, t) = e −at cosh(bkt) + a bk sinh(bkt) , (C4) where we have used the symmetries of the hyperbolic functions.At times t 0 much larger than a characteristic time of the system (t 0 1/bk), we can assume sinh(bkt 0 ) = cosh(bkt 0 ).3c, but with the effective long-time expression for ρ(k, t).The expansion shown in Eq.C5 is around t0 = 2/a, and except for the critically damped mode, we obtain a good agreement with the original form for large t.

FIG. 3 .
FIG.3.Fourier coefficients for density and longitudinal momentum in the underdamped and overdamped case, respectively.
FIG. 4. Simulation setup for molecular dynamics simulations of confined fluids.One lateral dimension is much larger than the gap height.

FIG. 5 .
FIG.5.Autocorrelation functions of the real part of density and momentum modes.The first column (a-d) corresponds to transverse momentum ( j ⊥ k), the second column (e-h) corresponds to longitudinal momentum ( j k), and the third column (i-l) corresponds to density correlations.Panels within the same row share the same wavelength.Solid lines show results from molecular dynamics (MD) calculations, and dashed lines are the corresponding fits.For simulations of the confined system, all wavelengths shown here are in the underdamped regime (470.6σ< λcrit).

1 τ 1 FIG. 6 .
FIG.6.Relaxation times obtained from fits to the autocorrelation functions of momentum and density fluctuation as a function of the wavelength for bulk and confined systems with α = 0.75.For all three subfigures, simulation results for the bulk system are shown as blue discs, and results for the confined system are shown as green diamonds.Panels a and b show the theoretical prediction for bulk and confined systems as dash-dotted and dashed lines, respectively.In panel b, the relaxation times according to the isothermal theory, Eq. (24b), are modified for short wavelengths by considering thermal effects, which is illustrated as a dotted line.Panel c shows thermal relaxation times obtained from density autocorrelation functions, which have quadratic wavelength scaling for both bulk and confined systems, as long as walls are rigid.Open diamonds in panel c show results from simulations with thermal walls, which indicate a similar transition as in a and b.

FIG. 7 .
FIG.7.Effective sound period and velocity of sound as obtained from longitudinal momentum correlations for bulk and confined systems from simulations with α = 1.5.In panel a, we observe an overall shift to longer sound periods for the confined fluids, which increases with wavelength, and panel b compares the sound velocities directly.For the confined system, the prediction based on the dispersion relation Eq.(20) adequately describes the deviation from the bulk.The critical wavelength marks the transition from underdamped to overdamped behavior at λcrit = 675 σ.
FIG.9.Time evolution of a small, Gaussian-shaped density perturbation in a slit geometry.The gap height in panel a is large enough to accommodate sound waves.Panel b shows the propagated distance ∆x over time t, rescaled by the critical wavenumber kcrit and the long-wavelength sound attenuation time τ = lim k→k crit τ , respectively.Therefore, we recover the sound speed cT = 1/τ kcrit from the slope of the curve.In panel c, the gap height is one order of magnitude smaller, which leads to a purely diffusive behavior.Panel d shows the variance σ 2 over time, rescaled by the squared critical wavenumber and sound attenuation time, respectively.Hence, the slope is proportional to the corresponding diffusion coefficient Dρ = h 2 c 2 T /12ν.Only the r.h.s. of the symmetric profiles is shown in a and b.

FIG. 10 .
FIG.10.Shear attenuation time for varying wall-fluid interaction and comparison to the bulk (a).The dash-dotted blue line illustrates a fit to the bulk relaxation time, and dashed lines show the wavelength-dependent shear relaxation time Eq.(24a) fitted to the data from confined systems.The long wavelength limit of the relaxation time determines the effective gap height, and since the geometric gap height is known, we can directly compute the slip length b using the definition of κ. Results of this fit procedure are shown in (b) for varying fluid wall interactions α, and are compared to non-equilibrium calculations to obtain the slip length.The dashed lines in (b) are a guide to the eyes.

FIG. 11 .
FIG.11.Same as Fig.3c, but with the effective long-time expression for ρ(k, t).The expansion shown in Eq.C5 is around t0 = 2/a, and except for the critically damped mode, we obtain a good agreement with the original form for large t.

TABLE I .
Theoretical expressions for the characteristic time scales for wave propagation and decay in bulk and confined systems.

TABLE II .
Material parameters obtained for a supercritical (T = 2.0 /kB, ρ = 0.452 σ −3 ) Lennard-Jones (LJ) fluid obtained from molecular dynamics (MD) simulations, and from the NIST database for supercritical Argon.Derived quantities such as e. g. the sound attenuation coefficient Γ may depend on values from both sources.All reported parameters are used for the predictions in Figs.6 and 7.