Molecular Mechanics of Liquid and Gas Slip Flow

By taking into account the inhomogeneity introduced by the presence of a solid boundary, slip-flow theory extends the range of applicability of the venerable Navier–Stokes description to smaller scales and into the regime where confinement starts to be important. Due to the inherently atomistic nature of solid–fluid interactions at their interface,slip flow can be described, at least in principle, predictively at this level. This review aims to summarize our current understanding of slip flow at the atomistic level in dilute gases and dense liquids. The discussion extends over the similarities and differences between slip in gases and liquids, characterization and measurement of slip by molecular simulation methods,models for predicting slip,and open questions requiring further investigation.


INTRODUCTION
Slip at a solid-fluid interface is said to occur when the tangential velocity of the fluid is not equal to the velocity of the solid it is in contact with (Lauga et al. 2007).It is typically a result of the inhomogeneity introduced by the presence of the solid and is typically captured by replacing the traditional no-slip boundary condition (Lauga et al. 2007) with a slip boundary condition.Despite its relative simplicity and certainly because of it, slip flow is an invaluable theoretical tool.It extends the traditional Navier-Stokes-based hydrodynamic models into the regime where confinement starts to become important while requiring only a change in the boundary condition, thus enabling researchers to continue using a well-understood model for which sophisticated solution methods, be they analytical or numerical, already exist.
By relaxing the requirement of no slip at the boundary, slip alleviates a kinematic restriction to flow, generally leading to flow rates higher than the nominal amounts expected in the absence of slip.In some cases, this flow enhancement can be significant; for example, in the limiting case of pressure-driven flow in nanotubes, it can completely dominate the response to the driving force (Kavokine et al. 2021), resulting in significantly larger flow rates or, conversely, requiring a significantly lower driving force for the same flow rate ( Joseph & Aluru 2008, Falk et al. 2010).In other words, fundamental understanding of slip-flow mechanisms is interesting from both a scientific and a practical point of view.
This review focuses on the fundamentals of slip flow, taking in particular an atomistic point of view.With this in mind, we consider simple fluids in the absence of electrostatic and non-Newtonian effects.A broader discussion of slip can be found in the review by Lauga et al. (2007), in which additional physical phenomena leading to or affecting slip, such as wall roughness and apparent slip due to electrokinetic phenomena or dissolved gas, are also discussed.
Transport in liquids and gases is fundamentally different.In dilute gases transport is essentially kinetic; particle interactions can be neglected, except during collisions (Vincenti & Kruger 1965, Cercignani 1988, Bird 1994).As a result, the internal fluid scale for hydrodynamic purposes is the molecular mean free path, ℓ.Deviations from the no-slip Navier-Stokes description are expected when the Knudsen number, Kn = ℓ/L, where L denotes the system characteristic scale, is no longer negligibly small.Typical values of ℓ for gases of practical interest at atmospheric pressure lie in the range of 50 to 100 nm, making slip corrections necessary for L ≲ 10 µm.The accuracy of the slip-flow description increasingly degrades as Kn increases, with its limit of applicability typically taken to be Kn ≈ 0.1, although this clearly depends on the amount of error that can be tolerated.For smaller length scales, a kinetic approach is required, which typically consists of solving the Boltzmann equation (Bird 1994, Homolle & Hadjiconstantinou 2007, Sone 2007, Radtke et al. 2011, Sadr & Hadjiconstantinou 2023) in the transition regime 0.1 ≲ Kn ≲ 10 or by approximating transport as ballistic for Kn ≳ 10.
Transport in dense liquids is dominated by intermolecular interactions.As a result, the internal scale for hydrodynamic purposes is the atomic/molecular size itself, σ , which for a simple liquid is of the order of 0.3 nm.Accordingly, deviations from the no-slip Navier-Stokes description are expected for L ≲ 10 nm.On the other hand, as recently elaborated on by Kavokine et al. (2021), for L ≲ 1 nm the continuum description is being confronted by structuring (Gravelle et al. 2013, Wang & Hadjiconstantinou 2017) discreteness and nonlocal effects.Slip flow performs a very important function in bridging the gap between these two length scales, since it eliminates the need for the development of another, intermediate, transport description.
As in the case of constitutive relations, despite fundamental differences in transport mechanisms, slip boundary conditions in gases and liquids are remarkably similar, at least to leading order.For isothermal, unidirectional flow of a Newtonian fluid in the x-direction parallel to Slip length, β, measured from the physical wall location, z = 0, shown here in the liquid case; the latter is defined by the locus of solid atom centers (gray) in direct contact with the liquid atoms (blue).In the interest of simplicity, the wall is shown as a single layer of atoms.The hydrodynamic wall location, z w , is discussed in Section 3.1.2.The distance between the first liquid layer and the wall, z min , is discussed in Section 3.2.1.Figure adapted with permission from Wang & Hadjiconstantinou (2017); copyright 2017 American Physical Society.
a stationary, flat solid boundary (wall) at z = 0 with normal in the positive z-direction (see Figure 1), the slip boundary condition is typically written in the form 1.
The slip parameter β can be interpreted as the distance inside the boundary at which the extrapolated velocity profile reaches the wall velocity, and for this reason it is also termed the slip length.
The canonical case of no slip is obtained in the limit β → 0. Note that the above relation is an already simplified form of the more general relation originally proposed by Navier (1823), where λ denotes the solid-fluid friction coefficient, xz = −µ(∂u x /∂z) is the x-z component of the stress tensor in the fluid, and µ is the fluid viscosity.Equation 2 is consistent with that obtained using general thermodynamic arguments (Bedeaux et al. 1976), which also provide for the presence of a thermal creep component.Thermal creep is briefly discussed in Section 2.2.2.Equation 1 directly follows from Equation 2 for β = µ/λ.
Recent careful experiments by Cross et al. (2018) demonstrate that by keeping the contribution of interface friction and viscosity separate, Equation 2 is more suitable for problems involving variable viscosity.
Significant research effort has focused on developing reliable methods for determining the value of the slip length (or friction coefficient) for various solid-liquid systems of interest.Despite considerable progress, several open questions remain, both in terms of robust and predictive methods for the calculation of these quantities and in terms of putting liquid slip flow on a firmer theoretical footing and on par with the dilute gas case.These topics are discussed in Section 3, following a presentation of the kinetic theory treatment of the dilute gas case in Section 2. Section 4 concludes with a summary and outlook for future work.

Background
Provided they are sufficiently dilute (Bird 1994), gases can be described by the Boltzmann transport equation (Vincenti & Kruger 1965, Cercignani 1988, Bird 1994), which governs the evolution of the single particle distribution function f = f (r i , c i , t), where r i = (x, y, z) denotes the position vector in physical space, c i = (c x , c y , c z ) denotes the molecular velocity vector, and t is time; here, we have assumed that no body forces are acting on the gas.Air at atmospheric pressure narrowly meets the dilute gas criteria (Bird 1994).
In the above, [f/t] coll denotes the collision operator, which describes the rate of change of the distribution function due to collisions in the gas.Within this description, intermolecular interactions are captured by the scattering law associated with binary collisions.The most well-known collision model is that of hard spheres (HS) (Vincenti & Kruger 1965, Bird 1994, Sone 2007).Due to the complexity associated with the HS collision operator, it is often replaced with simple models such as the Bhatnagar-Gross-Krook (BGK) (Bhatnagar et al. 1954, Cercignani 1988, Sone 2007) or the ellipsoidal-statistical model (Hollway 1966).An alternative that better represents the temperature dependence of transport coefficients of real gases has also been developed and is known as the variable HS model (Bird 1994).A survey of molecular interaction models and their resulting transport properties can be found in the review by Sharipov & Seleznev (1998).In the interest of simplicity, we discuss slip-flow theory in the context of the HS gas, which strikes a reasonable balance between fidelity and complexity.The HS gas mean free path is given by ℓ = ( √ 2ρπd 2 ) −1 , where ρ denotes number density and d the HS diameter.
The slip expressions discussed below are independent of the collision model, whereas the slip coefficients (nondimensional slip lengths, such as β/ℓ) are weakly dependent on the latter.On the other hand, slip coefficients are strongly dependent on the gas-solid interaction model.Within the Boltzmann equation framework, the latter is treated by a boundary condition on the distribution function at the physical wall location r w i .The most well-known such boundary condition is by Maxwell (1879), which assumes that the gas particle impinging on the boundary is diffusely reflected with probability α or specularly reflected with probability 1 − α.Diffuse reflection, also termed full accommodation, implies that the particle is reemitted from a Maxwell-Boltzmann distribution at the temperature and velocity of the boundary, T w and u w i , respectively, where Tw = 2k B T w /m, k B is the Boltzmann constant and m is the fluid (gas/liquid) atomic mass.The coefficient ρ g w can be thought of as the number density of an equilibrium dilute gas with the same temperature and velocity as the boundary; its value is determined by mass-flux balance at the boundary where dV c denotes the differential volume element in molecular velocity space.Using Equation 4to evaluate the first integral yields As one would expect, slip at the boundary increases as the degree of accommodation, α, decreases below 1.On the other hand, experiments suggest that engineering surfaces are most closely represented by fully accommodating models (Wenski et al. 1998, Arkilic et al. 2001, Sharipov 2003, Graur et al. 2009); for this reason, most results on slip coefficients can be found in this limit.Although some additional gas-surface interaction models have been developed (e.g., the Cercignani-Lampis model; Cercignani & Lampis 1971, Cercignani 1988), they have received significantly less attention, and as a result, they are less well characterized.In the interest of simplicity, in what follows we limit our discussion to the Maxwell boundary condition given by Equation 4.
The emergence of the Navier-Stokes description from the Boltzmann equation can be demonstrated using the well-known Chapman-Enskog expansion (Vincenti & Kruger 1965, Cercignani 1988).This expansion considers an infinite expanse of gas and thus provides no information about the gas behavior at the boundaries.Slip boundary conditions arise from the inhomogeneity introduced by the boundary presence.They can be rigorously derived using an asymptotic matching procedure (Sone 2007) for Kn â 1 in which the system is sufficiently large that a Navier-Stokes description can be found in the bulk of the physical domain and the boundary presence can be introduced as a perturbation to that solution.
Since the present interest in slipping flows is related primarily to microscale and nanoscale applications, in the interest of simplicity, we limit the discussion to low-speed, steady flows governed by the steady Stokes equation.Extensions to unsteady flows as well as those governed by the Navier-Stokes system of equations are briefly discussed once the theory has been presented.The thermal counterpart to the slip condition, the temperature jump, is obtained using an entirely analogous procedure (Sone 2007).Jump boundary conditions for phonon-mediated heat transfer in solids have also been derived using such a procedure (Peraud & Hadjiconstantinou 2016).
According to the asymptotic analysis, away from the boundary (in practice, a distance of the order of a few mean free paths) the Boltzmann solution is fully described by the Stokes system of equations.Close to the boundary, however, the distribution function associated with the Stokes flow field cannot satisfy the kinetic boundary condition, introducing additional kinetic effects and an inner, Knudsen layer component to the solution.As a result, the solution of the linearized Boltzmann equation for the gas hydrodynamic velocity, is written as a superposition of a Stokes component and a Knudsen layer correction, with u KN i → 0 away from the boundary.Slip boundary conditions emerge from the matching procedure involving the inner and outer solutions at each order of the expansions and

First-Order Slip
In the present section we summarize the salient results of such an analysis to first order in the Knudsen number.For more details the reader is referred to the monograph by Sone (2007).

Slip relations.
The asymptotic analysis sketched above yields the boundary conditions and for the components u S,0 i and u S,1 i , with the (nondimensional) Stokes set of equations governing the dynamics at each level and the wall velocity expanded in the form (Sone 2007, Aoki 2001).In Equation 12, β = β/ℓ, θ = θ/ℓ represent nondimensional slip lengths, or slip coefficients; moreover, and where T is the temperature field, which can be obtained by the steady heat conduction equation subject to analogous jump boundary conditions (for details, see Sone 2007).Additionally, ri = r i /L, ti denotes the unit vector tangential to the boundary, ni denotes the unit vector normal to the boundary, and c r = √ 2k B T r /m is the most probable molecular speed based on the reference temperature T r .Here we note that the no-slip boundary condition emerges as the boundary condition of the zeroth-order description, with slip appearing as a correction at higher orders.
It is possible to verify that the above set of governing equations and boundary conditions is equivalent to O(Kn) to the Stokes set 15.
with the "implicit" slip boundary condition as usually presented in the fluid dynamics literature.Given that the value of the viscosity is not variable (see discussion below), Equation 16is consistent with Equation 2 and the thermodynamic result given by Bedeaux et al. (1976).

Discussion.
In analogy with the treatment of the kinetic boundary condition on the distribution function (e.g., Equation 4), within the slip-flow description, the boundary condition is applied at the physical wall location.Moreover, impermeability requires that u S,1 i | w ni = 0.That is, no slip normal to the wall exists to this order; note that u w i ni = 0 for steady flow (Aoki 2001).Crucially, within slip-flow theory, the constitutive relation of the gas, including the value of the viscosity, remains unchanged.At the same time, no Knudsen layer correction exists for the stress.In other words, even though the Stokes component of the solution is unable to match the Boltzmann solution for the velocity field in the vicinity of boundaries, it is able to correctly capture the stress field in the same region.Adjusting the gas viscosity to enforce a match in the velocity field in the wall vicinity is not justified on practical grounds (due to its possible corrupting effect on the stress field, which requires no correction) or on theoretical grounds: A discrepancy between the Stokes flow field and the Boltzmann solution close to the boundary is required by the theory; if closer agreement is required in this region, a Knudsen layer correction should be used (see below).
The final term in Equation 12 is known as thermal creep and it describes the tendency of the gas to slip along a solid surface in the presence of a temperature gradient in the same direction (from cold to hot).This can be seen by noting that the term G S,0 i ti can be written as (c r /T r )(∂T 0 w /∂ ri )t i , since at the zeroth order (in analogy to Equation 11) we have T S,0 = T 0 w at the boundary (Aoki 2001).
The structure of the slip-flow Equations 12 and 16 is independent of the gas molecular model (Sone 2007).On the other hand, coefficients β and θ are mildly dependent on the gas model, while β in particular is strongly dependent on the degree of surface accommodation.Due to historical reasons, they are typically approximated by expressions that date back to Maxwell's first attempt at estimating them (Maxwell 1879) [e.g., β(α) = (2 − α)/α].Although those expressions are perhaps surprisingly accurate as estimates, they are not quite exact.For example, for a HS gas in contact with a fully accommodating wall, β = 1.11 and θ = 0.573 (Sone 2007); for the whole range of accommodation, Sharipov & Seleznev (1998) report that β can be approximated reasonably well by the relation β and Wakabayashi et al. (1996) report that θ varies between 0.439 (α = 0) and 0.573 (α = 1).Differences between various gas models (or approximations like the one by Maxwell) become particularly important in interpreting experimental results, especially those aimed at inferring the accommodation coefficient from slip coefficient measurements.Since differences between models can be sufficiently large to affect the inferred value of the accommodation coefficient beyond typical experimental error magnitudes, such experiments should always make clear which expression connecting slip to the accommodation coefficient was assumed in this process (see, for example, the work by Graur et al. 2009).
More recent work has shown that using the more realistic Lennard-Jones (LJ) molecular model leads to essentially the same slip coefficients (Sharipov 2011).A thorough compilation of slip coefficient data and their dependence on the molecular model and the gas-surface interaction model and parameters can be found in the review by Sharipov (2011).
The asymptotic solution approach also provides the following result for the Knudsen layer correction (Sone 2007): where η is a stretched coordinate normal to the wall, i.e., η = χ/ℓ, where χ denotes the coordinate normal to the wall.In other words, remarkably, the Knudsen layer correction depends only on velocity gradients of the Stokes solution at the wall and can thus be superimposed to the Stokes solution once the latter has been calculated.As in the case of slip coefficients, functions h 0 (η) and h 1 (η) depend weakly on the molecular model.They have been calculated primarily for the BGK and HS gas and are typically available in tabular form.
Figure 2 illustrates these concepts for the case of Couette flow of a HS gas subject to diffuse reflection (α = 1) at the system boundaries (z = ±L/2).The latter are moving with velocities ±U in the x-direction, with L = 10ℓ.The Stokes solution with the slip boundary condition as given by Equation 16 yields u even though strictly nonzero at distances as large as ∼10ℓ, the 90% thickness of this layer, defined as the distance from the wall at which the Knudsen layer correction has decayed by 90% of its value at the wall, is ∼2ℓ.These features are essential to the success of the slip-flow description.
Figure 2c highlights the important fact that the gas constitutive relation remains unaltered close to the boundary ( xz = −2µU/L(1 + 2 βKn)).Overall, the slip-corrected Stokes description remains accurate in the slip-flow regime and, as we remark below, remarkably robust beyond this regime.

Second-Order Slip
In the present section we summarize the salient results of the second-order asymptotic analysis.For more details the reader is referred to the monograph by Sone (2007) and the more recent work by Takata & Hattori (2012a).

Slip relations.
Solution to the second order proceeds along the same lines: u S,2 i is governed by a Stokes set of equations subject to a corresponding set of boundary conditions.At this order, these boundary conditions include boundary curvature effects (Sone 2007, Takata & Hattori 2012a), which are suppressed here in the interest of simplicity.Specifically, for a flat solid wall, the boundary condition for the tangential component of the flow velocity is given by

Discussion.
As in the first-order case, close to the boundary, each second-order slip term (including the terms due to boundary curvature not shown above) is supplemented by a secondorder Knudsen layer correction that is scaled by the same field gradients.The associated kinetic boundary layers are available only in tabular form (Sone 2007).At this order, gradients normal to the boundary lead to a slip flow in the normal direction.Zero mass flux across the boundary is restored by the Knudsen layer correction (Sone 2007).
For the case of a HS gas, for a fully accommodating wall, δ = 0.71 and φ = 1.25 (Takata & Hattori 2012a).Information on the dependence of these coefficients on the molecular potential, accommodation coefficient, and scattering model is limited.
In the case of unsteady flow, additional terms appear in Equation 18as well as in the Stokes set of equations (Radtke et al. 2012, Takata & Hattori 2012a, Takata et al. 2012).A comprehensive exposition of the unsteady problem and a discussion of the main differences from the steady case can be found in Takata & Hattori (2012a).
As in the case of the first-order theory, equivalent results up to O(Kn 2 ) can be obtained by solution of the Stokes set of Equation 15 subject to the "implicit" second-order slip boundary condition written here for isothermal unidirectional flow.Beyond convenience, in some cases the selection of slip relation (e.g., Equations 11 and 12 versus Equation 16or Equations 11, 12, and 18 versus Equation 19) may come down to a choice of class of boundary conditions that are preferable; this is especially true in the case of the slip condition given by Equation 19, where the highest-order derivative is of the same order as that of the governing differential equation, raising questions about posedness.
When using second-order boundary conditions such as Equation 18or 19, we must recall that the contribution of first-order Knudsen layers may become important.One specific example of practical importance is the calculation of the flow rate across the domain (Hadjiconstantinou 2006).Taking pressure-driven flow in a 2D channel with fully accommodating walls at z = ±L/2 as an example, the flow rate from the Stokes component of the flow field is given by where dP/dx denotes the pressure gradient along the channel.On the other hand, when calculating the flow rate to O(Kn 2 ), we need to include the first-order Knudsen layer contributions (Cercignani 1964, Hadjiconstantinou 2003, Takata & Hattori 2012b), namely, where in taking the limit of the integral of the function h 0 to infinity, we assumed that ℓ â L.
By the argument just presented, the contribution of second-order Knudsen layers to the flow rate would be of the order of Kn 3 .The above expression for Q can be written in the form (Cercignani 1964, Hadjiconstantinou 2003) ) dz, 22. where ψ = ∫ ∞ 0 h 0 (η)dη = −0.168.For the pressure-driven flow example, this expression leads to in which an effective second-order slip coefficient ( δ + ψ) appears that is different from the actual second-order coefficient δ.This is particularly important for the interpretation of experiments aimed at determining the second-order slip coefficient via flow rate measurements.Specifically, measurement of the flow rate enhancement due to slip will return the value of the effective slip coefficient δ + ψ and not the value of the actual slip coefficient ψ.
Figure 3 shows a comparison between first-order slip flow (boundary condition given by Equation 16), second-order slip flow (boundary condition given by Equation 19), and direct Monte Carlo simulations of Poiseuille flow of a HS gas between accommodating walls for Kn = 0.05, 0.1, and 0.2.The comparison shows that both first-order and second-order theories remain accurate even at the nominal limits of their applicability (Kn ≈ 0.1).Even more remarkable is the robustness exhibited well beyond these limits; for example, the second-order flow rate prediction (Equation 23) remains accurate to approximately 2.5% at Kn = 0.2.

SLIP IN DENSE LIQUIDS
Compared with the dilute gas case, our understanding of liquid slip is significantly less advanced.This is due to the persistent and multibody intermolecular interactions, which make theoretical treatments significantly more complex.As a result, molecular dynamics (MD) simulation (Allen & Tildesley 1989, Frenkel & Smit 2002) has been used extensively to characterize slip.In what follows we review various approaches for calculating slip based on MD simulations, as well as a variety of molecular-mechanics-based approaches for modeling this phenomenon.In the interest of simplicity we revert our notation to that of Figure 1; namely, we consider a liquid flowing in the x-direction, parallel to a stationary flat wall with normal pointing in the z-direction.

Molecular Simulation
Broadly speaking, MD simulations can be used to measure slip via either an equilibrium or a nonequilibrium approach.Both approaches are discussed below.To aid the discussion, we introduce the generalized LJ potential (Alosious et al. 2019) which serves as a canonical interaction potential in fundamental studies of slip.Here, u ij (r) denotes the interaction energy between atoms i and j as a function of their distance r, while C ij modifies the degree of attraction between atom i and j from the standard LJ case (C ij = 1) (Allen & Tildesley 1989).In what follows, we use subscript s to denote solid (wall) atoms and their properties and l to denote liquid atoms and their properties.Unless otherwise stated, physical quantities are reported in nondimensional units using σ ll as a characteristic distance, ε ll as a characteristic energy scale, and τ LJ = √ mσ 2 ll /ε ll as the characteristic time (here m = m l ).More realistic empirical potentials for simulating common material interfaces of interest (e.g., water in carbon nanotubes; Perez-Hernandez & Schmidt 2013) have been developed and are currently in use.Unfortunately, the performance of such potentials remains subject to a trade-off between accuracy and computational cost, meaning that inexpensive predictive calculations of slip for general solid-liquid systems currently remain out of reach.For a brief discussion of limitations associated with MD modeling of fluids under nanoconfinement, the reader is referred to the review by Kavokine et al. (2021).

Nonequilibrium method.
Slip lengths in MD simulations can be measured in the presence of a velocity gradient by extrapolating the velocity profile to the location where the fluid velocity is equal to that of the wall, as originally defined (see Figure 1).This approach is termed the nonequilibrium method.MD simulation is uniquely suitable for slip characterization in the nonlinear regime where approaches based on linear response cannot be used.On the other hand, these simulations need to be performed with care, as viscous heating, associated thermostats, or other nonlinear effects may contaminate the result (Martini et al. 2008a, Pahlavan & Freund 2011).
As we discuss in more detail below, typical applications of interest fall within the linear regime, which can be accessed in nonequilibrium MD simulations using low driving forces, with the beneficial side effect of minimizing the abovementioned artifacts.Although statistical noise makes MD simulation at small driving forces expensive (Hadjiconstantinou et al. 2003), current computational resources make such calculations possible on modern workstations.
An additional phenomenon that complicates slip measurement in the nonequilibrium method is liquid layering at the solid-liquid interface (Wang & Hadjiconstantinou 2017).At this time no firm theoretical guidance exists on the effect of layering on the local velocity profile (in analogy to the Knudsen layer in the dilute gas case, for example).Using Couette flow simulation to measure the slip length has the advantage that the velocity gradient term in Equation 1 can be measured in the liquid bulk, thus excluding the layered region (Sendner et al. 2009, Hadjiconstantinou & Swisher 2022b) where local properties vary as a result of the liquid density variation (Wang & Hadjiconstantinou 2018).
A final source of uncertainty in such simulations is related to the wall location and is a result of the finite size of atomic interactions.Measurement of the slip length with respect to the physical wall location (z = 0 in Figure 1) yields a different result than measurement with respect to the first liquid layer (FLL) or the hydrodynamic wall location (z = z min and z = z w in Figure 1, respectively); the difference in the first case is of the order of one atomic diameter (see Section 3.2.1),while in the second case it may be more (see Section 3.1.2).

Equilibrium method and Green-Kubo formulations.
As a dissipative transport coefficient, the friction coefficient at the solid-liquid interface can also be measured with a Green-Kubo formulation (Frenkel & Smit 2002), which requires only equilibrium simulations.In the formulation originally proposed by Bocquet & Barrat (1994), leading to where we use the superscript GK to denote the result of the Green-Kubo calculation.In the above expressions, angled brackets denote ensemble average, A denotes the interface area, and F x (t) denotes the force exerted by the solid on the fluid in the x (slip)-direction.
A severe limitation to the accuracy of this approach arises from the difficulty in identifying the plateau in the integral ∫ ∞ 0 ⟨F x (t )F x (0)⟩dt in finite systems (Bocquet & Barrat 1994, de la Torre et al. 2019).This issue is typically sidestepped by approximating (Bocquet & Barrat 1994) the plateau with the first peak of Oga et al. (2019) has shed some light on the error committed by this approximation in finite systems.Using several assumptions, they show that can be approximated by where the parameters 0 , t 1 , and t 2 contain information about λ GK and the two relevant characteristic system timescales, namely, the GK and the viscoelastic relaxation timescales.The value of the friction coefficient (within the model assumptions) can be obtained from the expression after determining the parameters 0 , t 1 , and t 2 by fitting traces of (t) to MD data.Using this model, Oga et al. ( 2019) calculated the error associated with the first-peak approximation as a function of the ratio of the GK to the viscoelastic system timescales.In particular, they show that the error becomes small when this timescale ratio is large, with a ratio greater than approximately 100 yielding an error of less than 5%.For systems that do not meet the required ratio of timescales (as determined by the desired error tolerance), they recommend using their model for determining the slip length.
Figure 4 shows a typical simulation result for an LJ fluid.Using Equation 28, which provides a remarkably good fit to the MD data, shows that, in this case, using the peak value of ∫ t 0 ⟨F x (t ′ )F x (0)⟩dt ′ underestimates λ GK (and thus overestimates β GK ) by approximately 4% compared with the result from Equation 29.The discrepancy typically increases with increasing λ GK .
Running Green-Kubo integral, Molecular dynamics (MD) simulation result for ρ = 0.68, T = 1.5, ε sl = 0.6, and C sl = 0.6.Simulation parameters and setup as in Hadjiconstantinou & Swisher (2022b).Blue dots denote MD data, and the red line denotes the best fit of Equation 28 to the MD data (with t 1 = 3.52, t 2 = 0.0483, and 0 = 0.324).Approximating λ by the (first) peak of the MD trace, as shown by the dashed gray line, results in an estimate approximately 4% smaller than that obtained from Equation 29.

The hydrodynamic wall location.
In our discussion so far, the location at which the slip boundary condition is to be applied coincided with the physical wall location (z = 0 in Figure 1).GK theory, on the other hand, predicts that the effective solid-liquid interface at which the slip boundary condition is to be applied is located inside the liquid, some distance from the physical wall location.This distance, also termed the hydrodynamic wall location (Bocquet & Barrat 1994), is given by where the x-z component of the liquid stress tensor (at equilibrium) is given in terms of atomistic quantities by xz = ∑ i m i c xi c zi + ∑ f xi z i .In this expression, i runs through all liquid atoms, c ji denotes the velocity of atom i in direction j, m i denotes the mass of atom i, and f xi denotes the x-component of the force on liquid atom i, including the forces exerted by the solid onto the liquid.
Calculation of z GK w by Equation 30 is quite challenging.In addition to Equation 30 being significantly more sensitive to noise (Bocquet & Barrat 1994) than Equation 26, a reliable approach for resolving the plateau dilemma for this ratio of GK integrals has yet to be developed.As a result, the hydrodynamic wall location has remained largely uncharacterized, with studies utilizing the GK formulation typically failing to report this quantity.Unfortunately, as shown recently (Hadjiconstantinou & Swisher 2022b), neglecting z GK w cannot be theoretically justified for moderate and especially small slip lengths, where it can be of the same order as, if not larger than, β.

Hydrodynamic wall location via tangential force balance at the wall.
The limitations associated with evaluating Equation 30 can be avoided by using the tangential force balance at the wall,  Comparison of slip measurements from equilibrium and nonequilibrium molecular dynamics (MD) simulations under a variety of conditions.The Green-Kubo prediction, β GK , as well as associated uncertainty is denoted by the dark red line and light red shading, respectively; the slip measured at z = 0 by nonequilibrium MD simulations, β, is shown by gold dots, while blue dots denote the measurements at the hydrodynamic wall location z = z GK w (β + z GK w ).For more details about the MD simulations, see Hadjiconstantinou & Swisher (2022b).
to provide a different relation for evaluating z GK w .This approach (Hadjiconstantinou & Swisher 2022b) uses measurements of the solid-liquid force in nonequilibrium MD simulations to determine z GK w as the location at which u x (z = z GK w ) satisfies the above equation.This method was extensively validated by Hadjiconstantinou & Swisher (2022b).Figure 5 shows typical results from this work, where nonequilibrium slip measurements from Couette flow were compared with equilibrium measurements of slip obtained from Equation 26.Nonequilibrium measurements are shown in two forms: as measured at the physical wall location, denoted by β, and referred at the GK hydrodynamic wall location, given and denoted by β + z GK w .This shift is possible because the flow profile in Couette flow is linear.As can be seen from Figure 5, only β + z GK w is in agreement with β GK ; that is, only when the hydrodynamic wall location associated with β GK is taken into account can the equilibrium and nonequilibrium measurements be brought into agreement.Duque-Zumajo et al. (2019) developed an alternative GK approach that addresses (Camargo et al. 2019, de la Torre et al. 2019) the plateau problem discussed above.Following Hansen et al. (2011), their analysis is based on the equation governing the dynamics of a liquid slab of thickness of a few atomic diameters immediately next to and in contact with the solid.

An alternative Green-Kubo formulation.
The final result of their analysis is where and 35.
In the above, V is the liquid volume and of the liquid-only stress tensor (at equilibrium): that is, fxi denotes the x-direction component of the force on liquid atom i from all of the liquid atoms (but not the solid atoms).
With the exception of the upper limit of integration, the expressions for µ D and λ D are the standard expressions for the liquid viscosity and the solid-liquid friction coefficient; G D can be interpreted as a correction to the liquid viscosity in the wall vicinity due to the wall forces (P.Español, personal communication).Although Equation 32 appears to suffer from the same limitation as Equation 30, namely, requiring evaluation of the ratio of GK integrals, according to Duque-Zumajo et al. ( 2019), Equation 32overcomes the plateau problem because, by design, both the numerator and the denominator exhibit a polynomial long-time tail, which leads their ratio to reach a plateau at some intermediate timescale, denoted here by τ D , from which the slip length can be calculated more reliably.
It can be shown (Camargo et al. 2019) that the set of Equations 32-35 is equivalent to Equations 26 and 30.The difference between Equation 32 and the standard relation β GK = µ/λ GK can be understood by noting that β D is based on a modified hydrodynamic wall location z D w = z GK w − G D β GK /µ.This translation of the hydrodynamic wall location in exchange for an equivalent change in the slip length is of course possible because both theories assume that the velocity profile is linear (in the wall vicinity).We also note that z D w is still to be determined by a ratio of two GK integrals (Camargo et al. 2019), making it susceptible to the limitations described above.

Slip Models
In this section we review atomistic models of slip at the solid-liquid interface.We start the discussion with a review of liquid layering at the solid-liquid interface, a well-known yet largely uncharacterized phenomenon.

Layering at the liquid-solid interface.
As in the case of the molecular velocity distribution in a gas, the wall presence introduces an inhomogeneity in the structure of a dense liquid in the wall vicinity.This inhomogeneity manifests itself as a region that extends approximately 3-5 atomic diameters from the wall in which the liquid structure is layered; depending on the solid-liquid system and its thermodynamic condition (see below), density variations in this region can be very appreciable.
As the degree of confinement increases, the size of this interfacial region increases relative to the full system size.It can thus affect both the amount of liquid in the system (Wang & Hadjiconstantinou 2015) and the local liquid properties (Giannakopoulos et al. 2012, Wang & Hadjiconstantinou 2018), and perhaps most importantly, momentum and energy exchange with Typical layered density profiles for a liquid density ρ ave = 1.The blue line indicates Wa = 1.2, and the red line indicates Wa = 0.6; h FL denotes the width of the first liquid layer (see text).Adapted with permission from Wang & Hadjiconstantinou (2017); copyright 2017 American Physical Society.
Despite the importance of layering in confined systems, most research work has taken a phenomenological approach.Here we note some recent work focused on the development of continuum-based computational approaches (Raghunathan et al. 2007, Mashayak & Aluru 2012, Mashayak et al. 2015) for accurately describing fluid structure under nanoconfinement.
A typical density profile of a layered liquid as observed in an MD simulation is shown in Figure 6.Such profiles are in good qualitative agreement with experimental results obtained using X-ray scattering (Huisman et al. 1997, Cheng et al. 2001).The width of the depletion region in the immediate wall vicinity (denoted by z min in Figure 1) can be estimated quite accurately by noting that it is primarily a result of the very strong short-range repulsion between the solid and the liquid.MD simulations show that due to the very steep form of the repulsive potential, it can be reliably determined as the distance from the wall at which the liquid atom potential energy due to all wall atoms is equal to zero; for a wall made of a single layer of atoms (e.g., a graphene sheet), this leads to z min = (2/5) 1/6 σ sl (Wang & Hadjiconstantinou 2017), and for a 3D wall this leads to z min = (2/15) 1/6 σ sl (Sendner et al. 2009).Both of these results are given for the case C sl = 1.
An operational definition of the degree of layering based on the ratio of the first maximum and following minimum in the density profile was proposed by Wang & Hadjiconstantinou (2017).MD simulations were then used to show that the wall number Wa = ρ w σ 2 sl ε sl /k B T , where ρ w denotes the wall number density, is the nondimensional number controlling the degree of layering as it represents the balance between solid-liquid attraction and thermal energy.Significant layering is typically expected when Wa is no longer very small, although the exact value depends on the wall details and definition of significant layering.We also note that the definition for Wa given above is for the case of a wall comprising a single layer of atoms and thus σ 2 sl ρ w represents a nondimensional areal wall density; for multilayered walls, a properly nondimensionalized volumetric wall density should be used.
Further insight into the layering phenomenon was presented in Wang & Hadjiconstantinou (2017), which pursued a theoretical investigation of layering on the basis of the Nernst-Planck equation, governing the density distribution of the liquid as determined by a balance between interatomic interactions and thermal energy.Referring to the geometry of Figures 1 and 6, the Nernst-Planck equation for the system of interest can be brought (Wang & Hadjiconstantinou 2017)  where ρ = ρ(z)/ρ ave , U sl denotes the liquid particle potential energy due to interactions with the wall, U ll denotes the liquid particle potential energy due to interactions with all other liquid particles, and ϵ = ρ ave ε ll /ρ w ε sl .In the above, ρ ave represents the average liquid density in the system (number of particles divided by the volume accessible to the liquid; also approximately equal to the liquid density away from the interface for large systems).
When scaled as shown above, Equation 36 can be used to motivate the role of Wa in controlling the degree of layering, with ϵ affecting the finer features of layered-liquid density profiles.Asymptotic analysis in the limit of small ϵ can be used (Wang & Hadjiconstantinou 2017) to show that typical quantities of interest, such as the areal density in the FLL, FL ; the fluid density in the FLL, ρ FL ࣕ FL /h FL ; and the width of the FLL, h FL , can be expressed in polynomial powers of ρ ave .
One surprising and at the same time convenient result is that the enhancement factor C = FL /(h FL ρ ave ) (excess number of particles in the FLL compared with the absence of layering) is to a large extent independent of ρ ave and a very weak function of temperature for ρ ave ≳ 0.7.In other words, over a wide range of temperatures and densities 1.5 ≲ C ≲ 2.5, a surprisingly small variation given the range of behaviors observed for different system parameters.A qualitative explanation follows by noting that the sharp peaks in the density profile are typically accompanied by lower minimum values of the density, leading to some cancellation (for an example, see Figure 6).
Much remains to be understood about the layering phenomenon and its effects on transport at the nanoscale in both simple and complex liquids.

Molecular-kinetic theory of liquid slip.
A simple, approximate model that synthesizes the main elements governing liquid slip on a molecularly smooth solid surface can be obtained (Wang & Hadjiconstantinou 2019) by modeling slip as a thermally activated process and using an extension of the Eyring theory of reaction rates (Eyring 1936, Hirschfelder et al. 1955) to model it.Evidence for the validity of the thermally activated process approximation is provided by the success of the related molecular-kinetic theory (MKT) of Blake & Haynes (1969) applied to the motion of contact lines, the direct validation of the model discussed below (see Figure 7), and the careful work by Lichter et al. (2004) and Martini et al. (2008a), who showed that slip is a thermally activated process at sufficiently low shear rate.Blake & Haynes's model in particular is widely accepted in the contact-line literature (Bonn et al. 2009) having been found to be in good agreement with experimental data (for a discussion, see Blake 2006) and MD simulations of contact-line motion (Blake et al. 1997, de Ruijter et al. 1999, De Coninck & Blake 2008, Wang et al. 2019).
This model considers the drift of liquid atoms in the FLL in contact with the solid.As such, it does not incorporate any notion of a hydrodynamic wall location, thus effectively assuming that the hydrodynamic wall location is the FLL.Although the distance between the FLL and the physical wall location is less than one atomic diameter, in the interest of completeness we distinguish between quantities measured at the two locations by referring to the slip velocity measured at the first layer of liquid atoms as u FL x and denoting the associated slip length β FL .According to the MKT, u FL x can be written (Brochard-Wyart & de Gennes 1992, Wang & Hadjiconstantinou 2019) as the difference between particle forward and backward hopping rates over the energetic barrier landscape they sample in the wall vicinity, namely, 37.
This leads to where τ 0 is the jump timescale, V is the potential barrier associated with the jump of length l, and ω is the modification to the potential barrier associated with the jump due to the presence of shear.
In writing Equation 38it is assumed that the potential barrier landscape V can be decomposed into a term capturing solid-liquid (V sl ) and a term capturing liquid-liquid (V ll ) interactions.Finally, by writing ω = −1 FL xz l = µ −1 FL γ l, where γ denotes the shear rate, we obtain the following expression for the slip velocity: Figure 7 shows that a scaled form of this expression, namely, u FL x = u 0 sinh( γ / γ0 ), can capture MD simulation results for slip as a function of shear rate from a variety of studies.Moreover, the values of coefficients u 0 and γ0 required to fit the MD data (Wang & Hadjiconstantinou 2019) implied reasonable values for the parameters τ 0 exp (V/k B T) and l.The comparison includes the data of Thompson & Troian (1997) Comparison between Equation 39, denoted by the gray dashed line, and molecular dynamics (MD) simulation results denoted by symbols [green symbols denote all five sets of data discussed by Thompson & Troian (1997), purple symbols denote the data of Martini et al. (2008a), and blue symbols denote the data of Wang & Hadjiconstantinou (2019)].Adapted with permission from Wang & Hadjiconstantinou (2019); copyright 2019 American Physical Society.Annu. Rev. Fluid Mech. 2024.56:435-461.Downloaded from www.annualreviews.orgAccess provided by 173.76.102.191 on 01/23/24.See copyright for approved use.
to explain the highly nonlinear dependence of the slip velocity on the shear rate.In the above, β 0 denotes the low-shear-rate slip length and γc denotes a critical shear rate.It immediately follows from Equation 39 that in the small-shear-rate limit γ ≪ 2 FL k B T /(lµ), Equation 39 linearizes to Equation 1, with Here we note that the natural appearance of the shear stress as the driving force makes this approach physically consistent with the general form of Equation 2. We also note that the low-shear (linear regime) condition can also be written as γ ≪ ρk B T /µ by noting that FL ∼ ρh FL ∼ ρl.Although the presentation above has followed the work of Wang & Hadjiconstantinou (2019), Blake and colleagues (Blake 1993, Blake et al. 2015) have independently discussed extension of his contact-line model to the general slip case and arrived at similar results.Related ideas can also be found in the work of Tolstoi (1952), Ruckenstein &Rajora (1983), andBowles et al. (2011), while the earliest discussion on the use of rate processes to describe slip at an interface involving a fluid appears to be the one found in the book by Glasstone et al. (1941).
Blake proposes eliminating the term involving the liquid-liquid interaction potential barrier using µ ∝ v −1 l exp(V ll /k B T ), where v l is the liquid flow volume (Blake 2006), and identifying the solid-liquid interaction with the work of adhesion V sl ∝ γ sl (1 + cos θ ), where γ sl is the liquid surface tension (in contact with vapor) and θ is the solid-liquid contact angle.The strong dependence of slip on the work of adhesion for the contact-line case has been extensively discussed by Blake in several publications (see, e.g., Blake 2006) and more recently used in calculations that include both contact-line and Navier slip (Keeler et al. 2022).A discussion of the earlier work on this topic by Tolstoi can be found in Blake (1990).
Although elegant in its simplicity, the MKT introduces several approximations, especially in its application of transition state theory (TST) (Hänggi et al. 1990) and its use of several unknown proportionality coefficients as well as atomistic-level lumped parameters, such as the flow volume and the jump length, which make predictive calculations for real materials of practical interest challenging.For a discussion of the connection between Eyring's reaction rate theory and TST in the context of slip (albeit for contact lines), the reader is referred to Blake & De Coninck (2011).

Detailed mechanistic models: Frenkel-Kontorova and Fokker-Planck descriptions.
Here we describe work aiming to construct a more detailed mechanistic model of slip and thus address some of the criticisms associated with the MKT discussed in the previous section.Such a model can be developed by writing the equation of motion for liquid atoms in the FLL in contact with the solid and subject to a uniform shear rate γ .Referring to the geometry of Figure 1, the equation describing motion parallel to the wall for a liquid atom in the FLL can be written in the form where overdots denote differentiation with respect to time, η sl is an atomistic-level solid-liquid friction coefficient (Martini et al. 2008b), A is the area per atom at the solid-liquid interface, f l→ x denotes the force exerted by all other liquid atoms on the atom of interest (beyond the shear force µ γ A), and f s→ x denotes the force resulting from the periodic energetic barrier associated with the layered solid.A reasonable approximation for the latter is f s→ x = − f 0 sin(2π x/a), representing the largest term in a Fourier decomposition of the liquid-solid interaction potential, with a denoting the lattice spacing in the solid (Steele 1973, Barrat & Bocquet 1999, Lichter et al. 2004, Martini et al. 2008b).
On the basis of their MD simulation results, Lichter and collaborators (Lichter et al. 2004, Martini et al. 2008b) constructed a detailed model for f l→ x , which includes nearest-neighbor interactions between liquid atoms in the FLL and the fluctuations in the number of atoms in this layer.Using a linear spring approximation of the form f l→ x = k △ x for those forces, where △ denotes the discrete Laplacian operator, transforms Equation 41 into a Frenkel-Kontorova (FK) equation (Braun & Kivshar 1998).Accounting for the motion of liquid atoms normal to the boundary (in and out of the FLL) results in a variable-density FK equation.Unfortunately, this equation is quite complex and could only be solved by simulation.Simulation results demonstrate qualitative agreement with in-house MD data (Lichter et al. 2004, Martini et al. 2008b).
The observation by Lichter and collaborators that at low shear rates, slip appears to be a thermally activated process suggests that a simpler approach might be reasonable, especially if it enables numerical or even analytical solutions.In one such approach (Hadjiconstantinou 2021), the term f l→ x is modeled as a Wiener process (Risken 1989).Using this simplification, Equation 41 can be written in the following rescaled form, with where ma, and d n = τ 2 f 0 /2πma.In the above, δ(t) is the Dirac delta function, τ = ντ LJ is a characteristic timescale, and ν is a constant of order 1.The single-particle distribution function n = n (x n , v n , t n ) corresponding to the dynamics described by Equations 42 and 43 is governed by the Fokker-Planck equation (Risken 1989) where v n = ẋn , describing Brownian motion in a periodic potential.Macroscopic slip corresponds to steady solutions of this equation on the interval 0 ≤ x n < 2π and can be calculated from the particle mobility Numerical and analytical results (Hadjiconstantinou 2021) show that the slip predicted captures all the salient features observed in MD simulations, namely, a linear regime of constant slip length at small driving forces, a rapid rise in slip for moderate/large forcing, and a plateau in the slip length for sufficiently high forcing, as originally discussed by Martini et al. (2008a).
The ability of this model to accurately capture MD simulation data was evaluated by Hadjiconstantinou (2021).Considering the dependence of slip on shear rate shows that a variety of MD data, including those from the work of Thompson & Troian (1997), can be fitted by solutions of Equation 44 with accuracy comparable to that achieved in Figure 7; the reader is referred to Hadjiconstantinou (2021) for details on the fitting process.Of particular interest was the observation that the MD data are best fitted using ζ n â 1 ( n ∼ 1), for which Equation 44 lends itself to asymptotic solution.Specifically, several semianalytical results are possible (Hadjiconstantinou 2021) in the limit ζ n / √ n ≪ 1; here we give a particularly simple result for the slip length in the limit of small forcing ( f n /ζ n â 1) and large potential barrier (d n / n → ∞): Although strictly asymptotic, this result is in agreement with the exponential dependence of slip on the normalized barrier height also seen in Equation 40arising from the simpler MKT treatment.

Models arising from partial evaluation of Green-Kubo expressions.
In this section we explore how the GK formulation of Section 3.1.2can be used to provide a model that connects slip to the underlying physical processes as well as system quantities that are accessible to experimental measurements.We begin by rewriting Equation 26 in the form where The importance of this rearrangement is that it reduces the contribution of the GK integral into two distinct factors, namely the mean square wall force (MSWF) ⟨F 2 x ⟩ and the timescale I β , which can be analyzed and potentially modeled separately.MD simulations show (Falk et al. 2010, Hadjiconstantinou & Swisher 2022b) that ⟨F 2 x ⟩ is strongly dependent on the solid-liquid interaction and fluid state, while the timescale I β is a very weak function of both.In what follows we discuss these two main ingredients in isolation.

The timescale I β .
In early work (Barrat & Bocquet 1999, Priezjev & Troian 2006), timescale I β was argued to be proportional to the diffusive relaxation of the in-plane densitydensity correlations at the liquid-solid interface.This result followed from an assumption that the density-density correlation function can be decoupled into a product of a diffusive relaxation and a static correlation function (Barrat & Bocquet 1999).This approximate approach required parameterization of this diffusive relaxation using MD simulations, making it no more accessible than the original Equation 26.
More recent (Hadjiconstantinou & Swisher 2022b), as well as related (Hadjiconstantinou & Swisher 2022a), work has suggested the following alternative approximation I β = τ µ /D β , where D β is a constant and τ µ is the homogeneous liquid viscous relaxation timescale defined by where is the liquid high-frequency shear modulus (Heyes et al. 2019).This yields which takes a step toward making the dynamics within the GK integral more transparent.
Figure 8 shows a comparison between Equation 51 and the MD results of Hadjiconstantinou & Swisher (2022b), with D β = 1.05.This value was chosen as the one that gives best overall agreement between Equation 51 and the MD results; G ∞ was calculated using the analytical expression for an LJ liquid given by Zwanzig & Mountain (1965), and A −1 ⟨F 2 x ⟩ was taken from MD simulations.The agreement observed in the figure is very encouraging given the variety of conditions and solid-liquid interaction parameters explored in the validation (Hadjiconstantinou & Swisher 2022b).In the above comparison, liquid properties in Equation 51 were evaluated at the bulk density (ρ ave ) instead of at the local density in the FLL (ρ FL ).This can be perhaps explained by the observation that the two densities are directly related, as discussed in Section 3.2.1.Validation of Equation 51 under a variety of conditions for a Lennard-Jones liquid.Molecular dynamics (MD) simulation results and associated uncertainty (for β GK ) are depicted with blue dots, and the prediction of Equation 51 and associated uncertainty are depicted with the dark red line and light-red shading, respectively.For more details about the MD simulations, the reader is referred to Hadjiconstantinou & Swisher (2022b).

The mean square wall force ⟨F 2
x ⟩.The MSWF is an equilibrium quantity.Its strong dependence on the liquid-solid interaction and fluid state leads to these quantities having a large effect on the slip length.Simple scaling analysis suggests that ⟨F 2 x ⟩/A ∼ FL ε 2 sl .This scaling behavior has been largely validated by various MD simulation studies (Priezjev & Troian 2006, Sendner et al. 2009, Hadjiconstantinou & Swisher 2022b).Particular attention has been paid to the scaling β ∼ ε −2 sl (Priezjev et al. 2005), including the connection to the wetting contact angle via the relation ε sl ∼ (1 + cos θ ) (Sendner et al. 2009).We also note that β ∼ ε −2 sl is in qualitative agreement with β ∼ exp(−ε sl ) predicted by the MKT of Section 3.2.2 for moderate values of ε sl .
A more complete model for the MSWF was developed by Barrat & Bocquet (1999) by rewriting ⟨F 2 x ⟩ in terms of the density-density correlation function (Boon & Yip 1991).By approximating the solid-liquid interaction force in the slip direction in terms of the dominant term in a Fourier decomposition, they reach the expression where S denotes the structure factor (Boon & Yip 1991) of the liquid inside the FLL.This expression correlates slip data quite well (within a multiplicative constant) in a variety of settings (Priezjev et al. 2005, Priezjev & Troian 2006, Falk et al. 2010), although in some cases empirical corrections are necessary (Barrat & Bocquet 1999, Priezjev 2007).Although unquestionably useful for providing insight into the physics of slip flow, Equation 52 requires knowledge of quantities that require MD simulations from which ⟨F 2 x ⟩ or even β could have been measured directly.In other words, a model that connects ⟨F 2 x ⟩ to material properties easily accessible to experiments is still needed.
We close by emphasizing that Equation 51 approximates the GK estimate of the slip, which is associated with the hydrodynamic wall location z GK w .

DISCUSSION
One of the lessons learned from more than a quarter century of research in microfluidics and nanofluidics is that the (Navier-)Stokes description is remarkably robust, at least in this application domain.A signature of this robustness is the low rate at which the accuracy of this description degrades when used beyond its nominal limits of applicability (see, e.g., Hadjiconstantinou 2005Hadjiconstantinou , 2006;;Falk et al. 2010).Slip boundary conditions are instrumental contributors to the wide range of scales over which this description can be applied.Both dilute-gas and dense-liquid slip-flow theories provide some opportunity for misuse.In the dilute gas case, the presence of Knudsen layers means that (Navier-)Stokes theory alone cannot capture the flow field in the vicinity of the walls.To use the theory correctly, one has to resist the temptation to force agreement in this region using ad hoc constructs with adjustable parameters.In the dense liquid case, complexity is introduced by the concept of hydrodynamic wall location, which depends on the liquid state in a nontrivial manner and can be at a distance from the physical wall location that can be a significant fraction of the slip length, when the latter is small.Unless the hydrodynamic wall location is properly accounted for, slip measurements from different methods will not be consistent.
The mathematical sophistication of the dilute-gas slip-flow theory stands in stark contrast to the simplicity of the interaction model on which it was primarily built, which dates back to the first papers written on the subject.The sensitivity of slip on the gas-surface interaction, as evidenced, for example, by the significant increase in slip as accommodation decreases, only serves to highlight this deficiency.In other words, more realistic gas-surface interaction models need to be both developed and incorporated within the asymptotic theory described in Section 2.
In addition to contributing to our fundamental understanding of the atomistic processes governing slip in liquids, tractable models are desirable as substitutes for expensive MD calculations but also as a means of relating slip to quantities accessible to experiments, potentially sidestepping MD limitations associated with the fidelity of interatomic interaction potentials.The models discussed in Section 3.2 are useful first steps toward a complete theory of slip for dense liquids.It is also encouraging that some aspects of these models appear to translate to slip of more complex fluids (Niavarani & Priezjev 2008, Priezjev 2009, Hatzikiriakos 2012).
Most theoretical results discussed in this review have been validated by molecular simulations.This can be attributed, at least in part, to the challenges associated with experimentally accessing the relevant length scales for quantifying slip, related phenomena, and their origins, especially in the dense liquid case.Although the contribution of molecular simulation methods in guiding and validating model development is invaluable and a testament to their ability to probe atomistic length scales, ultimately experimental validation of theoretical or simulation results is necessary, as additional effects, not included in MD simulation or equivalent models, may be important (for example, see the conclusions in Secchi et al. 2016).
As discussed by Lauga et al. (2007), experimental investigations of slip can be performed by direct or indirect measurement, with the latter approach relying on a measurement of the effect of slip on some other quantity.A discussion of some of the experimental setups used to quantify slip can be found in the reviews by Lauga et al. (2007) and Bocquet & Charlaix (2010).These discussions show that, in general, indirect methods are more common, perhaps as a result of the flexibility they afford.On the other hand, converting the experimental measurement into a prediction for slip requires a carefully constructed, accurate model for the effect of slip on the measured quantity (Cottin-Bizonne et al. 2008, Secchi et al. 2016, Cross et al. 2018).Second-order slip coefficient measurement in dilute gases serves as a cautionary tale of the importance of the model connecting slip to the quantity measured in the experiment; as discussed in Section 2.3, measurements of the flow rate enhancement compared with a no-slip baseline provide information on an effective second-order coefficient that is different from the actual slip coefficient, a fact that has been neglected in a large number of published studies.
The GK-based model of Section 3.2.4 can be improved by finding a description for the MSWF in terms of quantities more easily accessible to experiments.The other main ingredient in this model, namely, the approximation I β = τ µ /D β , would benefit from further characterization of its limits of applicability and from methodologies for calculating the value of D β from first principles.
By taking a higher-order gradient into account, second-order slip provides a measurable improvement in accuracy in the case of dilute gases.Development of a rigorous second-order slip model for the dense liquid case may prove to be equally beneficial.

DISCLOSURE STATEMENT
The author is not aware of any biases that might be perceived as affecting the objectivity of this review.

Bubble Plumes in
Figure 1

Figure 2
Figure 2Couette flow of a hard sphere gas between fully accommodating walls.Comparison between slip-flow theory (gray lines) and direct Monte Carlo simulation (DSMC) (orange dots).(a) Flow profile; the two solutions differ due to Knudsen layers close to both walls.Inset highlights difference in the vicinity of z = 0.5L.(b) Knudsen layer in the vicinity of z = −0.5L:difference between the DSMC result and the slip solution (gold dots) compared with the prediction of Equation17, with h 0 (η) taken from Sone (2007) (blue line).Despite the long tail, its effective (90%) thickness is ∼2ℓ.(c) Shear stress; gas constitutive relation is unmodified by the wall presence and no Knudsen layer exists for the stress field.
Figure 5 Figure 6 Figure7shows that a scaled form of this expression, namely, u FL x = u 0 sinh( γ / γ0 ), can capture MD simulation results for slip as a function of shear rate from a variety of studies.Moreover, the values of coefficients u 0 and γ0 required to fit the MD data(Wang & Hadjiconstantinou 2019)   implied reasonable values for the parameters τ 0 exp (V/k B T) and l.The comparison includes the data ofThompson & Troian (1997), who proposed a shear-dependent slip length of the form Figure 8 Nature Silvana S.S. Cardoso and Julyan H.E. Cartwright p p Consequences for Personal Exposure Rajesh K. Bhagat, Stuart B. Dalziel, M.S. Davies Wykes, and P.F.Linden p p p p p p p p p p p p p p 405 Molecular Mechanics of Liquid and Gas Slip Flow Nicolas G. Hadjiconstantinou p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p 435 Multiscale Velocity Gradients in Turbulence Perry L. Johnson and Michael Wilczek p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p 463 Fluid-Elastic Interactions Near Contact at Low Reynolds Number Bhargav Rallabandi p pp p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p 491Learning Nonlinear Reduced Models from Data with Operator Inference Boris Kramer, Benjamin Peherstorfer, and Karen E. Willcox p p p p p p p p p p p p p p p p p p p p p p p p p p p p p 521 Flow Mechanics in Ablative Thermal Protection Systems Nagi N. Mansour, Francesco Panerai, Jean Lachaud, and Thierry Magin p p p p p p p p p p p p p p p 549 Fluid Dynamics of Airtanker Firefighting Dominique Legendre p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p of corrections to Annual Review of Fluid Mechanics articles may be found at http://www.annualreviews.org/errata/fluid 31.