Hydrodynamic slip in nanoconfined flows: a review of experimental, computational, and theoretical progress

Abstract

Nanofluidics has made significant impacts and advancements in various fields, including ultrafiltration, water desalination, biomedical applications, and energy conversion. These advancements are driven by the distinct behavior of fluids at the nanoscale, where the solid-fluid interaction characteristic lengthscale is in the same order of magnitude as the flow conduits. A key challenge in nanofluidics is understanding hydrodynamic slip, a phenomenon in which fluids flow past solid boundaries with a non-zero surface velocity, deviating from the classical no-slip boundary condition. This review consolidates experimental, computational, and theoretical efforts to elucidate the mechanisms behind hydrodynamic slip in nanoconfined flows. Key experimental methods, such as the surface force apparatus, atomic force microscopy, and micro-particle image velocimetry are evaluated alongside emerging techniques like suspended microchannel resonators, dynamic quartz crystal microbalance, and hybrid graphene/silica nanochannels, which have advanced hydrodynamic slip characterization at the nanoscale. In addition to direct slip measurement techniques, methods like sum frequency generation spectroscopy, X-ray reflectometry, and ellipsometry are discussed for their roles in probing solid-liquid interfacial interactions, shedding light on the origins of hydrodynamic slip. The review also highlights the contributions of molecular dynamics simulations, including both non-equilibrium (NEMD) and equilibrium (EMD) approaches, in modeling interfacial phenomena and slip behavior. Additionally, it explores the influence of factors such as surface wettability, shear rate, and confinement on slip, emphasizing the interaction between liquid structuring and solid–liquid interactions. Advancements made so far have uncovered more complexities in nanoconfined flows which have not been considered in the past, inviting more investigation to fully understand and control fluid behavior at the molecular level.

---

Abdul Aziz Shuvo

Abdul Aziz Shuvo is a Ph.D. candidate in Mechanical Engineering at The Pennsylvania State University and an Assistant Professor (on leave) at Bangladesh University of Engineering and Technology (BUET). His research focuses on interfacial phenomena, electronics cooling, turbulence, and multiscale modeling. He earned both his Bachelor's (2019) and Master's (2022) degrees in Mechanical Engineering from BUET. Before his appointment at BUET, he briefly served as a lecturer at Ahsanullah University of Science and Technology, Dhaka (AUST). Later, he was promoted to Assistant Professor at BUET in 2022.

Luis E. Paniagua-Guerra

Dr. Luis E. Paniagua-Guerra is a Postdoctoral Scholar in the Mechanical Engineering Department at The Pennsylvania State University. He received his PhD and master's degree in mechanical engineering from The Pennsylvania State University and a bachelor's in mechanical engineering from the University of Guanajuato in Mexico. His major interests include the fundamental study of nanoscale thermal, hydraulic, and mechanical properties via atomistic level modeling, the study of solid–liquid interfacial phenomena, and the thermal characterization of liquid cooling solutions for electronics.

Juseok Choi

Juseok Choi is a Ph.D. candidate in the Department of Chemical Engineering at Pennsylvania State University, under the supervision of Prof. Seong H. Kim. He has studied the mesoscale structure of cellulose microfibrils in plants using sum-frequency generation (SFG) vibrational spectroscopy, with support from the Center for Lignocellulose Structure and Formation (CLSF). His current research focuses on the effect of the interfacial water structure at the solid–liquid interface on hydrodynamic slip behavior through experimental approaches.

Seong H. Kim

Prof. Seong H. Kim joined the faculty of Chemical Engineering at Penn State in 2001 after completing a Ph.D. in Chemistry at Northwestern University and conducting post-doctoral research at the University of California, Berkeley. He is currently a Distinguished Professor of Chemical Engineering and has courtesy appointments with the Department of Chemistry and the Department of Materials Science and Engineering at Penn State. Dr Kim's research interests lie in surface science and materials characterization. He has published more than 360 papers, and the h-index of his published work is 71. He recently finished writing a textbook on “Surface and Interface Analysis – Principles and Applications”, which is currently under production through Wiley and will come out in March 2025.

Bladimir Ramos-Alvarado

Dr. Bladimir Ramos-Alvarado is an Associate Professor of Mechanical Engineering and the principal investigator of the Interfacial Phenomena Lab (IPHEL) at the Pennsylvania State University (Penn State), he is a member of the Materials Computation Center at Penn State, and an Associate of the Institute for Computational and Data Sciences (ICDS) at Penn State. He also serves as a Subject Matter Editor for Applied Thermal Engineering (Elsevier). Dr. Ramos obtained a Bachelor's and a Master's degree, both in Mechanical Engineering, from the University of Guanajuato, Mexico, in 2009 and 2011, respectively. Later on, he continued his education at the Georgia Institute of Technology (Georgia Tech) where he obtained a PhD and a Master's degree in Mechanical Engineering in 2016. After brief Postdoctoral and Instructor appointments at Georgia Tech, Dr. Ramos joined the Department of Mechanical Engineering at Penn State University as an Assistant Professor in May 2017.

---

1 Introduction

Investigations on nanoconfined fluids date to the first half of the last century. Following early investigations of fluid transport in nanoconfinements,¹ the study of fluid flow at the nanoscale emerged as a new field of science called nanofluidics, a term coined in the 1990s.^2,3 Furthermore, the integration of nanofabrication with nanofluidic technology offers potential breakthroughs in the manipulation and control of fluids at the molecular level, which could revolutionize our approach to a variety of scientific and industrial challenges. This surge in interest is driven by the capability of nanofluidic devices to operate with high efficiency and sensitivity due to their minute scale and the unique properties of fluids confined to such small dimensions.

The application of nanofluidic devices could impact promising technologies, including biomedical applications, ultrafiltration and desalination processes, and ionic transport for energy conversion. In the biomedical field, nanofluidics may reduce the amount of genomic material and time for analysis of DNA.^4–6 Nanochannel devices have therapeutic applications on high-precision drug delivery systems;^7–10 while, nanoengineered fluidic devices enable low-cost cell analysis and disease diagnostics.^11,12 Similarly, the integration of nanoconduits in lab-on-a-chip systems allows for single-cell analysis, increasing the reliability of portable point-of-care medical diagnostic systems.^13,14 Beyond biomedical applications, nanoconfined fluid flows offer precise manipulation of ionic concentrations and enhance electrochemical-mechanical energy conversion in batteries.¹⁵ Such technologies also facilitate groundbreaking electrical energy production from salinity gradients.^16,17 Furthermore, the active control of ion charge concentration in nanoflows is applied to the development of fluid-based devices analogous to micro-electronics, such as nanofluidic transistors^18–20 and diodes.^21–24 Another main research venue of nanofluidics takes advantage of the filtration capabilities of nanotubes and nanopores.²⁵ Capillary devices with molecular-scale precision are used in ultrafiltration,^26,27 while nanometer-scale porous membranes are promising alternatives for seawater desalination.^28–32 In addition to the highly specialized applications discussed so far, the behavior of thin, confined fluids is essential in numerous industrial operations, mainly involving lubrication,³³ nanoencapsulation,³⁴ and nanofabrication.³⁵ Electrospinning technologies leverage the interaction between electrostatic forces and working fluids at the spinneret nozzle to create innovative nanofabrication methods to produce hollow, core–shell, and multichannel nanofibers, which are crucial for high-performance materials used in catalysis, drug delivery, and energy storage.³⁵

Along with the enhanced capabilities of nanofluidic devices comes a major increase in flow friction, which is an inevitable challenge. Fluid flow through minuscule confinements experiences a vast increase in flow resistance, per the classical hydrodynamics theory. The significance of friction in nanoconfined flows is attributed to the interfacial interactions as the surface-to-volume ratio increases at the nanoscale.³⁶ However, due to the dominant effects of interfacial properties such as wettability and surface roughness, fluid flow in nanochannels can be tailored to overcome the hydraulic limitations imposed by high confinement levels. For instance, experimental investigations of flow in carbon nanotubes (CNT) have reported flow enhancements up to five times higher than expected from the conventional continuum theory.^37–39 CNTs are atomically smooth surfaces with hydrophobic properties, a combination that hypothetically offers a reduced resistance to flow.

Nanoconfined flows have been described by classical fluid dynamics combined with atomistic modeling to account for interfacial interactions in nanochannels with diameters down to 1.4 nm.⁴⁰ The applicability of the continuum approach facilitates the creation of a theoretical framework to model fluid flow. However, the boundary condition compatible with the continuum approach remains a challenge, i.e., accounting for and quantifying hydrodynamic slip.^41–43 The earliest efforts to simulate flow through nanochannels resulted in flow velocities between 2 to 10 orders of magnitude higher than experimental measurements,⁴⁴ indicating an overestimation of the hydrodynamic slip. Sophisticated simulations, considering complex wettability interactions⁴¹ and long averaged simulated times⁴⁴ (∼100 μs) have been conducted in recent years, obtaining flow regimes comparable with the flow velocities reported in the experimental literature. However, a complete understanding of the mechanisms and significant parameters in nanoconfined flows is still needed.

In this literature review, the basics of hydrodynamics in nanoconfined flows are introduced in Section 2. A review of the experimental efforts aimed at explaining the different variables affecting hydrodynamic slip is presented in Section 3; the most relevant experimental techniques and the tendency towards higher resolution measurements are highlighted as well as the controversy regarding the effects of wettability and liquid structuring metrics on the hydrodynamic slip. Subsequently, a more comprehensive review of the modeling efforts for investigating the complex phenomena of nanoconfined flows is reported in Section 4. The different non-equilibrium and equilibrium methods for modeling fluid friction in nanochannels are revised and compared. Lastly, in Section 5 a more fundamental review of the underlying physics behind slip as well as the effects of wettability, shear rate, and confinement is presented. Fig. 1 illustrates a knowledge map outlining the contents of this review.

Fig. 1 Knowledge map of the literature review on nanoscale hydrodynamic slip.

2 Friction and hydrodynamic slippage in nanoconfined liquids

In macroscale fluid dynamics, the no-slip boundary condition has allowed experimental results to align with numerical and analytical models across various applications. Despite its widespread use, the no-slip boundary condition is phenomenological and not derived from fundamental physical principles.⁴⁵ In Navier's early work,⁴⁶ an alternative boundary condition that permits slippage has been suggested:where u_S and L_S represent the slip velocity and slip length, respectively. Here, z is the coordinate normal to the interface where the velocity gradient is assessed, and L_S is defined as the distance at which the linearly extrapolated velocity reaches zero. L_S = 0 denotes the no-slip condition. However, this slip boundary condition, as described in eqn (1), is also empirical and lacks a solid theoretical foundation.

Conversely, the governing equations of continuum fluid mechanics have a solid theoretical foundation. In the continuum framework, Newton's second law is applied to infinitesimal volume elements large enough to preserve the bulk values of the thermophysical and transport properties of the fluid. After applying the Newtonian corollaries to the relation between stresses and deformations, the famous Navier–Stokes (NS) equations arise. However, these equations break down when the flow system is reduced to dimensions comparable to the molecular size due to the uncertainty of the continuum assumption.⁴⁰

In the hydrodynamic regime where the NS equations are applicable, there is a distinct separation in length and time scales between bulk properties and surface effects.³ In bulk systems with particle quantities in the order of Avogadro's number, the molecular degrees of freedom can be described with only a few variables such as pressure, velocity field, temperature, etc., while the complexity of the transport phenomena is lumped into transport coefficients. Bocquet and Charlaix³ calculated the applicability range of the NS equations by determining the lower-scale limit for the concept of shear viscosity (η), which is represented by the Green–Kubo relation:

where V, k_B, T, 〈σ_xy(t)σ_xy(0)〉_eq are the system volume, Boltzmann constant, absolute temperature, and autocorrelation function of the off-diagonal component of the stress tensor, respectively. The validity of this equation assumes that the timescale of the stress-stress correlation function τ_σ is smaller than any hydrodynamic timescale. For instance, the relaxation time of momentum, τ_q = (υq²)⁻¹, where υ is the kinematic viscosity and q is a wave vector; thus, υq²τ_σ < 1 fixes the limit for timescales at confinements of size w larger than a viscous length scale, namely . For water τ_σ ∼ 10⁻¹² s and υ = 10⁻⁶ m² s⁻¹ at 20 °C, which yield a limit of w ≈ 1 nm, proving the robustness of the NS equations. Thus, for water, confinement levels in the nanometer scale can be modeled using the classical governing equations of fluid mechanics. Notably, Thomas and McGaughey⁴⁰ confirmed the bulk-like behavior of water flowing through CNT of diameters ∼1.4 nm.

3 Experimental investigations of hydrodynamic slip

Nanoconfined flowing liquids can be described by the equations of bulk hydrodynamics down to scales of approximately 1 nm (approximately three water molecular diameters). However, characterizing the boundary conditions for these flows remains challenging. The no-slip boundary condition has been successfully applied to macroscopic systems, where it remains phenomenologically valid. However, deviations from this condition arise at smaller scales where surface interactions dominate, casting doubt on its applicability at the nanoscale.

Accurately measuring interfacial properties in nanofluidics is challenging due to the dynamic interactions at liquid–solid interfaces and spatial resolution limitations.³ Recently, these experimental challenges have been addressed significantly. For instance, the surface force apparatus (SFA) and atomic force microscopy (AFM) can directly access the hydrodynamic forces with high sensitivity in a controlled environment. Additionally, micro-particle image velocimetry (μ-PIV) enables the visualization of the velocity field in real time, which is critical for slip measurements. This Section will discuss these primary methodologies—SFA, AFM, and μ-PIV—and their recent developments in measuring L_S in nanoconfinements. Moreover, Section 3.4 will explore emerging techniques for measuring L_S using suspended microchannel resonators (SMRs), dynamic auartz crystal microbalance (QCM-D), and hybrid graphene/silica nanochannel techniques. Along with L_S metrologies, surface characterization techniques, such as sum frequency generation (SFG) spectroscopy and X-ray reflectometry (XRR), will also be presented in Section 3.5. These techniques are useful in uncovering the origins of hydrodynamic slip behavior, particularly concerning the chemical interactions of interfacial water molecules.

3.1 Surface force apparatus (SFA)

The SFA measures the viscous force, F_h, between two surfaces submerged in a liquid with viscosity η, based on their separation distance h, and balanced by the cantilever restorative force, F_k,⁴⁷ as shown in Fig. 2(a). This viscous force is detected using piezoelectric materials, and the gap between the surfaces is determined through interferometry. Employing SFA in conjunction with multiple beam interferometry enables measurements with sub-nanometer accuracy. The expression for F_h is provided in eqn (3) where v_a represents the relative approach velocity of the surfaces, r is the radius of the sphere immersed in the liquid, and f* is a correction factor that compensates for hydrodynamic slip. When f* = 1, eqn (3) corresponds to the NS solution in the lubrication approximation. For surfaces of comparable characteristics, Vinogradova⁴⁸ derived the solution for f* as detailed in eqn (4).

Fig. 2 Schematic of the operational principles of (a) SFA, (b) AFM, and (c) μ-PIV. In (a) and (b), the radius of the spherical probe (r), surface separation distance (h), and probe approach velocity (v_a) are indicated as the parameters of the hydrodynamic drag force (F_h) and the restorative force of the cantilever (F_k), which are balanced for the hydrodynamic slip measurement.

Chan and Horn⁴⁹ investigated the drainage of three non-polar organic Newtonian liquid films between molecularly smooth mica surfaces using the SFA. Their results were in good agreement with the Reynolds theory of lubrication (no-slip boundary condition) for film thicknesses above 50 nm; however, for thinner films, an apparent enhancement of the liquid's viscosity was observed. Chan and Horn⁴⁹ indicated that as the film thickness decreases, the surface effects are more significant and produce a “solid-like” ordering in the liquid layers near the wall, causing an increase in the liquid resistance to shear. A modification to the formerly static SFA apparatus is reported in Luengo et al.⁵⁰ where a shear attachment was added to the original design, now allowing for dynamic rheological analyses. A series of transition regimes in the rheological behavior of polymer melts as well as a reduction of viscosity with film thickness were found using the modified SFA.

Baudry et al.⁵¹ used a similar version of the dynamic SFA to investigate the slippage of glycerol in contact with wetting and non-wetting surfaces. On a wetting cobalt surface, the no-slip boundary condition was observed; alternatively, when cobalt was coated with a thiol, the surface became non-wetting, and L_S ≈ 65 glycerol molecular diameters was measured. The liquid confined between wetting surfaces showed a constant viscosity as the confinement was varied; in contrast, for non-wetting cobalt-thiol, a reduction in the liquid viscosity was observed as the separation between surfaces decreased.

High shear rates have been observed to trigger bubble nucleation, leading to increased slip. Zhu and Granick⁵² observed large L_S with water and tetradecane on hydrophobic surfaces, noting that L_S escalated unboundedly as shear rates increased in their experiments, reaching up to the micrometer range. Furthermore, Cottin-Bizonne et al.⁵³ investigated the hydrodynamic boundary conditions of water and dodecane on both hydrophobic and hydrophilic surfaces, employing a dynamic SFA under shear rates up to 5 × 10³ s⁻¹. They found that the viscosity of water remained consistent with its bulk value in confinements as narrow as 10 nm. On hydrophilic Pyrex surfaces, neither water nor dodecane exhibited hydrodynamic slip; however, on hydrophobic surfaces, both fluids showed slip with L_S reaching approximately 20 nm.

The findings by Cottin-Bizonne et al.⁵³ present notable contradictions to those of Zhu and Granick.⁵² Cottin-Bizonne et al.⁵³ noted that the oscillation amplitude of pressure in their SFA experiments remained below the vapor pressure of the examined liquids, precluding cavitation, which might have been possible under the conditions used by Zhu and Granick.⁵² Furthermore, they hypothesized that discrepancies could arise from the contamination of surfaces by hydrophobic materials. Continuing this investigation, Cottin-Bizonne et al.⁵⁴ also explored potential experimental inaccuracies affecting their results. They highlighted that even small miscalculations in measuring the separation distance between surfaces could significantly impact the calculated L_S. Their further studies using water and water mixtures aimed to discern the impact of viscosity on different wettability surfaces, finding that surfaces with higher hydrophobicity exhibited greater L_S, though not exceeding 20 nm.

3.2 Atomic force microscopy (AFM)

AFM, as described in Fig. 2(b), employs the same approach and equations eqn (3) and (4) as the SFA. However, AFM probes smaller areas, which are determined by the size of the spherical bead attached to the AFM probe.³ Craig et al.⁵⁵ utilized an AFM setup to assess the drainage force in aqueous sucrose solutions between a spherical tip and a flat surface. They employed analytical models to determine variables such as L_S from the measured force–distance curves. They observed that the L_S varied with the fluid's viscosity and the rate at which the AFM cantilever approached (shear rate). At lower velocities, there was no slip, indicative of a “free system” behavior, whereas at higher velocities L_S ≤ 20 nm were reported. A notable limitation of AFM in measuring viscous forces in nanoconfined liquids is its sensitivity to deflections in the AFM cantilever caused by viscous drag.⁵⁶ To address this, Vinogradova et al.⁵⁶ developed multiple models aimed at curbing or even eliminating the increase in viscous drag as the speed of the AFM cantilever escalated. Further, Vinogradova et al.⁵⁷ engineered an AFM probe that reduced drag. Employing a data reduction method previously outlined by Vinogradova et al.,⁵⁶ they successfully identified the no-slip boundary condition and measured an L_S of 10 nm on hydrophobic surfaces.

Using AFM, Bonaccurso et al.⁵⁸ assessed the hydrodynamic force between hydrophilic mica and glass in the presence of aqueous solutions, measuring L_S of 8–9 nm on these surfaces, irrespective of the AFM probe approach speed. The experiments were conducted at high shear rates (10⁴ s⁻¹), which accounted for the observed hydrodynamic slippage on the hydrophilic surfaces. Honig and Ducker⁵⁹ explored the impact of rapidly approaching AFM probes on measuring viscous forces within sucrose solutions varying in viscosity. They observed no significant increase of L_S exceeding zero, even at shear rates up to 2.5 × 10⁵ s⁻¹, contrasting with the results from Bonaccurso et al.⁵⁸ when examining similar systems. However, their findings were in line with those of Vinogradova et al.⁵⁷ The discrepancies are thought to stem from the different methods used to measure the gap between the surface and the probe. Traditional AFM experiments calculate this distance using the combined displacement of the piezoelectric scanner movement and the cantilever deflection, while Honig and Ducker⁵⁹ derived the separation distance through the intensity of scattered evanescent waves. Bhushan et al.⁶⁰ investigated L_S for different surface conditions using the AFM tapping mode. The authors reported L_S values of 43 nm and 232 nm for hydrophobic and superhydrophobic surfaces, respectively. Maali et al.⁶¹ enhanced the design of commercially available AFM cantilevers by adapting them for improved acoustic excitation in liquid environments. They integrated an anti-reflective coated glass slide into the cantilever holder to reduce unwanted oscillation peaks. Subsequently, Maali et al.⁶² used this improved version of the AFM to measure the L_S of water on graphitic-carbon surfaces at a value of 8 nm.

In AFM experiments, instrumental uncertainties such as offsets and improper calibrations result in significant errors in determining L_S.^54,63 Such errors originate from the necessity of independently providing accurate values for the liquid viscosity, spherical AFM probe radius, and AFM cantilever spring constant when calculating L_S, as inaccuracies in these parameters can propagate and lead to incorrect slip estimations.⁵⁴ For example, Bonaccurso et al.⁵⁸ reported L_S of approximately 8 nm for a hydrophilic substrate (mica–water interface), while Maali et al.⁶⁴ and Zhang et al.⁶⁵ reported values close to zero. To overcome these discrepancies, a new data analysis model was proposed, assuming that the diameter of the microsphere is much larger than the slip length.^54,66,67

where h is the separation between the microsphere and the surface of the substrate, v is the approach speed of the sphere to the substrate, F_h is the hydrodynamic force acting on the tip, η is the viscosity of the liquid, and R is the diameter of the microsphere. Because ν/F_h is proportional to h + L_S, the slip length is determined by its intercept on the x-axis without knowing the viscosity of the liquid and the size of the microsphere. Ishida et al.⁶⁸ proposed an analytical alternative to eqn (5) adding the assumption that the L_S of the substrate is identical to that of the microsphere on the AFM tip. With their new analytical approach, they proposed that the L_S of a mica substrate with water is close to zero.

3.3 Micro-particle image velocimetry (μ-PIV)

The μ-PIV method, as illustrated in Fig. 2(c), involves tracking particle movement within a liquid flow constrained to microscale dimensions. This technique focuses image velocimetry analysis on areas close to the surface to sample the velocity profile and observe interfacial phenomena. Although theoretically applicable to nanoscale conduits, μ-PIV faces significant challenges in accurately tracking particles that are only a fraction of the size of the nanochannels. Additionally, as the resolution of μ-PIV improves, the measurements become noisier due to increased Brownian motion affecting smaller tracer particles.³ Consequently, μ-PIV is typically employed in microchannels, with efforts concentrated on enhancing resolution near the surface to detect the presence of hydrodynamic slip. Tretheway and Meinhart⁴⁵ utilized standard μ-PIV to investigate slip in hydrophilic and hydrophobic channels with a cross-section area A_{cross-section} = 30 × 300 μm². They used fluorescent polystyrene spheres with a diameter of 300 nm as tracers in the sampled region measuring 25 μm × 100 μm. They reported L_S = 1 μm for hydrophobic surfaces, while no-slip (L_S = 0) for hydrophilic conditions. These findings exceeded theoretical predictions but aligned with the expected effects of wettability.

Lumma et al.⁶⁹ enhanced the precision of velocity profiling within a 100 μm wide microchannel by using the fluorescence correlation spectroscopy (FCS) method to cross-correlate the fluorescence signals of clearly identified tracer particles with a diameter of 40 nm. This approach distinguishes between flow and diffusion effects, revealing L_S ranging from 0.2 to 1 μm. They recognized that their measurements of L_S might be larger, influenced by interactions between the liquid and surface and the repulsive forces among the tracer colloids. In another investigation using conventional μ-PIV, Ou and Rothstein⁷⁰ measured L_S of 7.5 μm in micro-ridges with ultra-hydrophobic surfaces and a shear-free configuration. Their results also matched independent experiments using a pressure drop calculation of slip. Joseph and Tabeling⁷¹ refined the μ-PIV method significantly, achieving near-simulation accuracy. They used fluorescent beads ranging from 100 to 200 nm in diameter within a microchannel measuring 10 μm × 100 μm. To precisely ascertain the wall position, they tracked the location of tracer particles that adhered to the walls, reporting L_S < 100 nm, with uncertainties comparable to the measured values.

3.4 Other techniques for direct slip measurement

In addition to the hydrodynamic slip characterization techniques described in Sections 3.1–3.3, there are methods to directly measure L_S. Collis et al.⁷² first proposed L_S measurements on individual gold nanoparticles immersed in water. They used suspended microchannel resonators (SMRs) introducing gold nanoparticles into “U-shaped” channels embedded in a cantilever. When a nanoparticle passes through the channel, it increases the inertial mass of the sensor. In the experiments, the flow at the particle's surface is closely related to the hydrodynamic boundary condition, and the slip at the particle surface is the result of the discrepancy between the measured and actual mass. Thus, L_S can be calculated by fitting the varying excitation frequencies of each vibrational mode against the corresponding mass discrepancies. The main advantage of using SMRs is the measurement without confinement conditions which can modify the nature of L_S. Unfortunately, further study of the hydrodynamic slip using SMRs has not been reported, but the authors suggested the measurement capability of particle wettability, particle crystal structure, particle surface functionalization, particle surface charge, system temperature, liquid viscosity, and polarity.

Xie et al.⁷³ developed a hybrid graphene/silica nanochannel for L_S measurements. The ratio of mass flow resistance between the silica and graphene sections was calculated using the meniscus movement in the graphene section and the corresponding capillary flow constants. The variation of mass flow resistance between graphene and the hybrid nanochannels is likely due to variations in slippage between the different sections of the hybrid nanochannel. This approach enables an indirect measurement of L_S in the graphene section of the nanochannel. In their result, the measured L_S = 16 nm is smaller than what has been estimated in the molecular dynamics simulation of graphene capillaries with pristine multilayered graphene (L_S = 60 nm).⁷⁴ They hypothesize that the observed variation of L_S is due to the functional groups and charges on the graphene surface during the deposition process.^75,76 This implies that there is growing interest in the relationship between hydrodynamic slip and interfacial chemistry.

Dynamic quartz crystal microbalance (QCM-D) is another emerging technique to study hydrodynamic slip. QCM-D measures changes in the resonant frequency of the crystal under oscillation, both before and after the deposition of mass onto the substrate; Zhdanov and Kasemo⁷⁷ reported the simulation-based observation of hydrodynamic slip on the surface of a QCM-D sensor. According to theoretical calculations, at low oscillation amplitudes, the amplitude of the substrate matched that of the central bead deposited on it, indicating a no-slip boundary condition. However, as the substrate oscillation amplitude increased, the oscillation amplitudes of the substrate and bead deviated from the theoretical prediction. Zhdanov and Kasemo⁷⁷ hypothesized that this mismatch is due to a transition from sticking to slipping at the interface between the substrate and the bead. Their additional finding that a lower transition amplitude occurred with a Ca²⁺-containing buffer (which changed the bead-support interaction) supports this hypothesis, suggesting that the frequency shift and energy dissipation in QCM-D may increase or decrease depending on the slippage between the substrate and the beads (or fluids). This approach could provide valuable insights into the influence of interfacial interactions on hydrodynamic slip, though further fundamental studies are needed to fully develop it.

3.5 Interfacial liquid characterization

The fundamental mechanisms of slip at solid–liquid interfaces are not fully understood, but theories and hypotheses suggest a significant influence of interfacial liquid properties, structure, and ordering on the nature of momentum transfer at solid–liquid interfaces. In this Section, analytical tools to probe the water/liquid interface at the molecular level will be reviewed. It should be noted that while the techniques presented here are not capable of directly measuring hydrodynamic slip, they can provide important properties linked to its origin and fundamental principles.

3.5.1 Sum frequency generation (SFG) vibrational spectroscopy. SFG vibrational spectroscopy is a non-linear optical process where two photons are combined at a surface and generate a new photon with its energy equal to the sum of two input photons. This process requires noncentrosymmetry. In the bulk liquid phase, all molecules are randomly moving, and such randomness is equivalent to centrosymmetry because the positive and negative directions on any axis are equivalent. Thus, bulk liquid cannot generate an SFG signal. In contrast, at solid–liquid interfaces, the randomness is broken, creating noncentrosymmetry; thus, SFG can detect interfacial molecular species without interference from the bulk phase molecules of the same species.

SFG has been extensively employed to investigate the interaction between water molecules and solid surfaces by analyzing OH stretching signals. One study on fused quartz in contact with liquid water revealed that all free OH groups at the silica surface are hydrogen-bonded to water molecules, which induces noncentrosymmetric ordering of water molecules near the surface. The distribution of disordered and ordered water structures in the interfacial region varies depending on pH.⁷⁸ The dipole direction of water at the solid–liquid interface flips by 180° when the pH of the aqueous solution crosses the isoelectric point of the surface.⁷⁹ This is significant as previous computational studies have demonstrated that the structuring of liquids at the solid–liquid interface plays a pivotal role in explaining the hydrodynamic slip behavior across different solid–liquid interfaces (see Section 5.2 for further information). Hence, the experimental characterization of liquid structure and orientation at the interface is crucial for bridging the knowledge gap between measurements and calculations of L_S.

When a solid surface interacts with a polar solvent such as water, it can be charged. This surface charge can arise through two mechanisms: either by the dissolution of surface groups into the contacting liquid and/or by the adsorption of ions from the solution. This process leads to the formation of electrical double layers (EDL), which can significantly influence solid–liquid friction. Using the SFG technique, Wei et al.⁸⁰ reported the relationship between pH, electrolyte concentration, and hydrogen bonding interactions at the interface. They observed a significant reduction of the SFG intensities of the OH stretching region (3000–3800 cm⁻¹) as increasing electrolyte concentration, and this drop was primarily attributed to the reduction of the number of water molecules orienting toward the solid surface in the EDL. Notably, numerical and theoretical models show that surface charges and salt concentration directly influence the hydrodynamic slip behavior. The reduction of L_S with surface charge was consistently reported in previous studies.^81–85 Geng et al.⁸³ highlighted that at high surface charge density, the L_S behavior is dominated by ionic interactions rather than solid–liquid binding forces. Rezaei et al.⁸⁵ conducted MD simulations to study the electro-osmotic flow of an aqueous NaCl solution on a charged silicon surface and observed a similar relationship between L_S and surface charge density.

Recently, Wang et al.⁸⁶ successfully detected SFG signals at graphene-water interfaces. Previously, isolating the graphene-water interaction was challenging due to interference from substrate-graphene signals because of the transparent nature of graphene. They employed a method of suspending graphene on the water surface, creating air-graphene-water interfaces under dry conditions to avoid signal interference from the air side. They observed that the OH peak appears at nearly the same frequency and amplitude as at the air-water interface, suggesting that graphene has only a weak effect on the organization of interfacial water. This result and approach open new opportunities to study the chemical interactions between graphene and water, which is of great interest for applications in microfluidic devices.

In addition to using spectroscopy metrologies for the characterization of solid–liquid interfaces, atomistic scale simulations are pertinent for investigating local transport properties, such as viscosity (η), the EDL thickness, and diffusion coefficients, particularly in nanoconfined electrolytes under varying electric fields and surface charge conditions.^87–89 For example, Masuduzzaman et al.⁸⁷ explored the effect of electric fields on nanoconfined aqueous electrolytes, showing that the EDL thickness increased while the local η decreased due to the bulk motion of counter-ions along the current flow direction. Further expanding this work, they also studied the impact of surface charge, revealing that the high local η in the first and second hydration layers results from strong electrostatic interactions and enhanced hydrogen bonding, which leads to a more ordered water structure.⁸⁸ This molecular ordering restricts mobility, increasing resistance to flow and thus increasing the local η compared to the bulk fluid (η_{first layer} > η_{second layer} > η_bulk). In addition, Masuduzzaman et al.⁸⁹ investigated the effect of the molecular interface position on the EDL thickness and hydrodynamic properties. Their findings showed that using the hydration layer as a boundary rather than the solid substrate's first atomic layer position resulted in better convergence toward the continuum assumptions. Furthermore, Ma et al.⁹⁰ conducted a comprehensive study on the relationship between water ordering and friction at the interfaces of water-TiO₂ and water-silicone using SFG and AFM. The wettability of TiO₂ and silicone substrates was systematically varied under different heating or plasma-treated conditions to assess its effect on L_S. Their results revealed a stark contrast: while hydrophobic TiO₂ substrates exhibited low friction, hydrophobic silicone surfaces demonstrated significantly higher friction. To investigate this potential discrepancy, the structuring of interfacial water molecules was analyzed through SFG. The spectra indicated that the interfacial water exhibited both loosely hydrogen-bonded “water-like” (with a peak at 3300–3600 cm⁻¹) and strongly hydrogen-bonded “ice-like” (with a peak at 3100–3300 cm⁻¹) structures, with the inhomogeneity in water structuring contributing to higher friction at the solid–liquid interface. Conversely, more uniformly ordered water structures (ice-like structures) in the first monolayer were reported to reduce friction on both TiO₂ (hydrophobic) and silicone (hydrophilic) surfaces. This suggests that the structuring of interfacial water, rather than wettability alone, plays a critical role in determining frictional behavior. To support these experimental findings, molecular dynamics simulations were conducted, showing that more ordered water structures reduced friction by decreasing hydrogen bonding and attractive interactions between the first monolayer and bulk water molecules. The reduction in hydrogen bonds and energy barriers facilitated a smoother interlayer movement, significantly lowering friction. In summary, Ma et al.⁹⁰ integrated both experimental and numerical approaches to provide deeper insights into the complex interplay between water structuring, friction, and hydrodynamic slip behavior.

3.5.2 X-ray reflectometry and ellipsometry. In Section 3.5.1, SFG spectroscopy was described as a surface-sensitive technique used to investigate the molecular water structure within the EDL region, which can influence the hydrodynamic slip behavior at solid–liquid interfaces. Studying the local liquid density at solid–liquid interfaces is also crucial to understand slip. Furthermore, computational investigations have demonstrated that the density depletion length can effectively explain hydrodynamic slip behavior; See section 5.2 for a detailed discussion. However, measuring the density profile of water molecules near the interface is challenging due to needing high surface sensitivity at the appropriate scale. In this section, X-ray reflectometry (XRR) and ellipsometry are introduced as unique optical techniques capable of measuring the density of water molecules near the surface.

XRR is used to obtain quantitative information on the electron density profile by observing changes in reflectivity (in-plane) of X-ray radiation at grazing incidence angles. In a multi-layer structure, reflected X-ray beams at each interface generate constructive or destructive interference patterns that can be compared to theoretical calculations based on the Fresnel equation. Using XRR, Mezger et al.⁹¹ determined the density profile of water/OTS(octadecyl-trichlorosilane)/SiO₂/Si layers, as shown in Fig. 3(a). They highlighted the depletion of water density at the interface with hydrophobic OTS layers, which was not detected by other experimental techniques. XRR was also used to investigate the surface of ionic liquid (1-methyl-3-octadecylimidazolium tris(perfluoroethyl) trifluorophosphate, [C₁₈mim]⁺ [FAP]⁻).⁹⁴ The coexistence of positively and negatively charged parts of the ionic liquid allows for structuring at the air–liquid surface. However, due to its structural complexity, the degree of molecular ordering by positively or negatively charged particles periodically oscillates with depth-decaying, showing the convergence of density away from the surface. This result demonstrates the capability of XRR to measure the density profile of ionic liquid multilayers at interfaces, as well as the limitations of general measurement methods for slip length (described in Sections 3.1–3.3), which typically involve fitting a single force–distance curve to calculate L_S.

Fig. 3 (a) Raw (red line) or step functioned (green line) density profile of bulk water/depletion/OTS/SiO₂/Si layers using XRR, (this figure has been adapted from ref. 91 from The National Academy of Sciences of the USA, copyright 2006), (b) schematic illustration of the experiment configuration of water/SiO₂/Nb₂O₅/BK-7 prism and refractive index profile used in the analysis of ellipsometry measurement (this figure has been adapted from ref. 92 with permission from Elsevier, copyright 2015), (c) density profiles obtained from MD simulations for different solid–liquid interactions (this figure has been adapted from ref. 93 with permission from American Chemical Society, copyright 2020).

Ellipsometry is also capable of detecting the density profiles of liquids near solid surfaces. It measures changes in polarization of the reflected light and relates them to the refractive index change as a function of distance. The refractive index profile can be related to the density profile via the Clausium-Mossotti relation. Wang et al.⁹² measured the refractive index profile of a water/SiO₂/Nb₂O₅/BK-7 prism under acidic and basic conditions. At pH 3, the surface is almost uncharged due to the isoelectric point of SiO₂, resulting in a negligible EDL, referred to as the absorbed layer in ref. 92 However, at pH 10, an additional layer representing the EDL is observed due to the more negatively charged silica surface by the deprotonation of the silanol groups, as shown in Fig. 3(b).

Although both XRR and ellipsometry lack chemical specificity, these techniques provide important structural properties (density and thickness) of liquid molecular layers near surfaces, which are critical for determining hydrodynamic slip. The experimental results from these techniques can be used to validate computational hydrodynamics simulations such as density profiles predicted from molecular dynamics as illustrated in Fig. 3(c).⁸⁸

3.6 Summary

In this Section, widely used experimental methodologies including, SFA, AFM, and μ-PIV are described for the direct measurement of L_S. SFA and AFM quantify the viscous drag force, which is balanced by the restorative force imposed on the surfaces compressing a liquid. μ-PIV can directly measure the velocity of liquid near a solid wall by tracking particle movement. Additionally, other approaches such as using SMRs, hybrid nanochannels, and QCM-D with varying flow resistances were presented. These techniques are capable of measuring L_S at the solid–liquid interface, but they also have limitations in providing detailed information on the fundamental mechanisms of hydrodynamic slip. Specifically, these dynamic measurement methods cannot provide any information on the liquid-surface interactions. Such information can be obtained with SFG, XRR, and ellipsometry. But, since these techniques work for the static condition, the link between the dynamic slip and the static interfacial structure remains elusive. The working principle and limitations of all these experimental techniques are summarized in Table 1.

Table 1 Summary of experimental techniques used to analyze the hydrodynamic slip behavior at the liquid–solid interface and the surface chemistry effect on liquid structure

Technique	Description	Limitation
Surface force apparatus (SFA)	• Measures the repulsive force as two surfaces approach one another	• Interpretation is model-based
	• Directly measures LS based on the rigid theories	• Requires highly controlled and smooth surfaces; limited to the chemistry control on mica
		• Limited to relatively small separations
Atomic force microscopy (AFM)	• Measures the repulsive force between a cantilever tip and a surface	• Interpretation is model-based
	• Sensitive to very small-scale interaction tip	• Sensitive to probe type and geometry
	• Directly probe LS on the surface	• Difficult to reproducibly control the chemistry of the probe surface
Micro-particle image velocimetry (μ-PIV)	• Tracks particle movement and velocity within a liquid flow confined to a microscale dimension channel.	• Lack of surface chemistry control
		• Limited resolution
		• Requires the use of an appropriate tracer
Microchannel resonators (SMRs), Hybrid nanochannel	• Measures mass changes of flows through a vibrating microchannel caused by the interaction at the surface boundaries (such as channel wall or nanoparticles)	• Lack of surface chemistry accountability
		• Requires complex fabrication and precise design
		• Complex interpretation
Dynamic QCM (QCM-D)	• Measures the frequency shift and energy dissipation of a vibrating quartz crystal as a liquid interacts with its surface	• Deconvoluting from the bulk fluid effect
	• Sensitive to very small changes in interfacial properties	• Insufficient studies on hydrodynamic slip
Sum frequency generation (SFG)	• Nonlinear optical technique by combining two laser beams to produce an SFG signal	• Does not provide a direct measurement of slip
	• Can provide molecular-level information about the alignment and behavior of liquid molecules at the solid–liquid interface	• Limited to surfaces and interfaces that are optically accessible
X-ray reflectometry (XRR)	• Measures the reflectivity of X-rays from a surface with high resolution and sensitive to atomic-scale structural changes	• Does not provide a direct measurement of slip
	• Can provide information about the density profile of a liquid near the solid interface	• Limited to well-defined, homogeneous surfaces
Ellipsometry	• Measures changes in polarization of the reflected light providing the density profiles of liquids near solid surfaces	• Does not provide a direct measurement of slip
		• Requires precise knowledge of the optical properties of the interface

---

Combining dynamic and static characterization methods as complementary probes is paramount to understanding the underlying mechanisms of hydrodynamic slip; however, such a combination of multiple experimental techniques for the same interfacial system has not been done systematically yet. In contrast, studies bridging the complex interplay between water structuring, friction, and hydrodynamic slip behaviors, have been done computationally, which is reviewed in Section 4.

4 Molecular dynamics (MD) simulations of nanoconfined flows

4.1 Non-equilibrium molecular dynamics (NEMD)

NEMD simulations have been applied to determine the friction coefficient and L_S for different solid–liquid interfaces. A key aspect of NEMD simulations is the application of a gradient across the computational domain to observe the system's linear response. The variables of interest are monitored over time, and once a quasi-steady state is achieved, transport or interfacial properties can be determined by fitting the data to a continuum-based model, typically involving a gradient. While the NEMD method is conceptually straightforward and closely mimics real experimental setups, it is significantly influenced by the size of the simulation box and the large gradients needed to extract statistically significant information, the latter due to the limited time scales accessible in MD simulations. Additionally, defining the location of the solid–liquid interface is paramount for performing the velocity extrapolation needed to compute L_S in the NEMD model. Karim et al.⁹⁵ reported that L_S and η calculations, based on a solid–liquid interface defined at the first hydration layer, closely matched predictions from the modified Hagen–Poiseuille equation.

For the analysis of nanochannel flow, NEMD simulations that resemble Couette and Poiseuille flows have been used. NEDM simulations consider flow between parallel plates, due to the restrictions imposed by the utilization of periodic boundary conditions in MD simulations. Generally, the simulations involve confining liquid molecules between solid walls and applying an external force to induce unidirectional motion of the liquid particles. To create Couette flow, one wall is moved at a constant speed while the other remains stationary,^96–101 or both walls are moved in opposite directions,^102–106 see Fig. 4. Couette flow is generally preferred for slippage and solid–liquid friction investigations due to its straightforward implementation and the constant shear rate observed over the whole liquid domain.

Fig. 4 NEMD methods for generating Couette flow, (a) moving walls in opposite directions and (b) moving one wall while keeping the other wall fixed.

Poiseuille flow simulations are used for more general purposes, such as studying the size effect on hydrodynamic slip¹⁰⁷ and the characterization of flow regimes in nanochannels.^108,109 Unlike Couette flow, the Poiseuille flow in NEMD can be modeled in several ways. One common approach involves applying an external force to the liquid particles within a defined region, often referred to as the inlet,^110,111 and allowing the system to reach equilibrium to achieve a parabolic velocity profile, as shown in Fig. 5(a). An optimization of this method proposed by Ge et al.,¹¹² considers the addition of an external force to the fluid particles in region B, and then the velocity is rescaled in region C to generate a constant inlet temperature every time the fluid moves between periodic boundaries from outlet to inlet, see Fig. 5(b). Another method to create a pressure-driven flow, introduced by Zhang et al.,¹¹³ involves confining fluid particles with three walls; one stationary at x = 0 and two moving walls advancing at a constant speed perpendicular to the stationary one, as illustrated in Fig. 5(c). This approach successfully produces a parabolic velocity profile and a linear correlation between the mean velocity and pressure drop. Lastly, Poiseuille-like flow can be generated by applying a constant force on every liquid atom in the desired direction of flow,^{107–109,114,115} see Fig. 5(d). Although this method does not resemble a pressure-driven flow but a body-force-driven flow, the outcome regarding the velocity profile and shear rates are similar.¹¹⁶

Fig. 5 NEMD methods for generating a Poiseuille flow: (a) inlet-driven flow where an external force and thermostat are applied in the same finite region (based on the methodology used in ref. 110 and 111), (b) optimized inlet-force-driven flow (based on the methodology used in Ge et al.¹¹²), (c) pressure-driven flow using three walls (based on the methodology used in Zhang et al.¹¹³), and (d) driving force applied to each liquid particle (based on the methodology used in ref. 107–109, 114 and 115).

To calculate L_S, a velocity profile is obtained from the particle trajectories. This is done by dividing the liquid region into segments (bins), each one recording the streamwise velocity of individual atoms, see Fig. 6. These velocities are then averaged per bin and over multiple timesteps to accurately delineate the velocity profile. The resulting profile is modeled based on the type of flow (shear-driven or pressure-driven) using either a linear or parabolic function. The slip velocity is determined by extrapolation to the wall or any other characteristic length (usually one molecular diameter from the wall). Finally, L_S is calculated as:

where Δu_s is slip velocity at the interface and ∂u/∂z is the velocity gradient calculated at the interface location z₀ where z corresponds to the direction perpendicular to the wall. If the observed velocity profile is plug-like, one can track the force and the slip velocity over several time steps and calculate the friction coefficient from λ = F_A/u_s, where F_A is the force acting on the solid walls per unit area and u_s is the slip velocity.¹¹⁷

Fig. 6 NEMD methodology to generate velocity profiles using discrete bins for (a) Couette flow (linear profile), and (b) Poiseuille flow (parabolic profile).

NEMD simulations of nanochannel flow are intuitive for implementation and analysis. Unfortunately, the high shear rates necessary to dampen the numerical noise, are significantly larger than the experimentally accessible shear rates. Consequently, the work done by the external forces produces viscous heating in the liquid particles, which is a collateral effect. To avoid this problem, many authors have applied the following methods: thermostating the liquid atoms while keeping the solid atoms frozen;^{96–98,113,118,119} thermostating all particles;^96,98,99,102 thermostating the fluid particles only in the perpendicular direction of flow (solid atoms are either mobile or inert);^{101,103,104,116,120–124} and thermostating only the solid atoms so they can act as natural heat sinks to remove viscous heating.^{107–110,114,115,125–129} Notably, most of the earlier NEMD nanochannel flow simulations employed the liquid thermostating method or all-particles-thermostating method due to the low computational cost of these strategies. Moreover, most of these investigations were concerned with the calculation of L_S in uniform temperature flows. Alternatively, the few early works where the solid thermostating method was applied were concerned with the relationship between heat dissipation and hydrodynamic slip.

The discussion of the thermostating approach is amply reported in the literature. Martini et al.¹⁰⁶ carried out multiple simulations to investigate the hydrodynamic behavior of a nanoconfined fluid under different shear rates. The main finding of this investigation was that when the liquid atoms are subjected to thermostating and the solid wall atoms remain frozen, L_S increases exponentially with shear rate; conversely, if the solid atoms are thermostated and allowed to vibrate, the L_S growth is bounded to a constant value at high shear rates. Ho et al.⁹⁹ indicated that either thermostating the liquid or the solid was consistent with experimental conditions of high and moderate heat dissipation, respectively. They also observed that the use of different thermostats on the fluid did not affect the velocity profiles.

In a more in-depth investigation, Bernardi et al.¹³⁰ indicated that thermostating could alter the flow physics if not applied correctly to a system where inhomogeneities exist, such as in nanochannel flow. Two scenarios were considered using a 2-D Couette flow: (1) thermostating the fluid particles while keeping the solid atoms "frozen" and (2) thermostating the solid atoms while allowing the liquid to heat up. The first finding was that the properties of the wall affect the density profiles. If the solid walls are allowed to vibrate, the liquid can slightly push the walls and make the channel effectively bigger. Moreover, thermostated fluids presented unrealistic shear distribution curves, which were not consistent with theoretical expectations. The effect of the wall atom dynamics was like that observed by Martini et al.,¹⁰⁶ showing that at low shear rates, the vibrating walls allowed larger slip than frozen walls. The wall dynamics effect was explained by looking at the elasticity of particle collisions using rigid and flexible walls (tethered by elastic springs). Depending on the collision angle, a rigid wall limits the direction of motion of a liquid particle post collision, whereas an elastic wall allows for liquid particles to easily move in the direction of flow after a collision (see Fig. 7(a) for the schematic of rigid and flexible wall types).

Fig. 7 (a) Schematic of rigid and flexible walls interacting with liquid particels in NEMD models, (b) temperature profile obtained by Shuvo et al.¹²⁹ (this figure has been adapted from ref. 129 with permission from AIP publishing, copyright 2024).

Yong and Zhang¹³¹ performed a series of MD simulations with three different thermostat setups: (i) applying the thermostat only to liquid atoms, (ii) applying it to both solid and liquid atoms, and (iii) applying it solely to solid atoms. They investigated the mechanical and thermal properties of the system for each setup. Their findings showed that using a thermostat on the solid walls resulted in a parabolic temperature profile, which aligned with the solution of the energy equation. However, deviations from theoretical expectations occurred when isothermal conditions were applied to the liquid atoms at high shear rates. Similarly, Shuvo et al.¹²⁹ applied thermostating to the walls only to mimic natural cooling, and observed a parabolic temperature distribution in the nanochannel consistent with theoretical expectations, as illustrated in Fig. 7(b). It can be concluded that allowing for wall heat dissipation, as it would naturally occur, is preferred over direct liquid thermostating.

As previously discussed, the most realistic manner of performing NEMD simulations of nanofluidics is by thermostating the wall atoms to allow the liquid to expel the excess viscous heating through the walls; however, a significant computational demand is inherent to this approach, i.e., the equations of motion must be solved for the solid particles too. Bernardi et al.¹³⁰ and recently De Luca et al.,¹³² proposed a novel methodology for conducting physically realistic simulations without solving the dynamics of the wall. This is a new thermostating method that uses virtual particles that only interact with the liquid while being tethered to their lattice site via an elastic spring model, see Fig. 8. The wall particles are rigid, while virtual particles serve as heat sinks. Although this method was proven to be efficient and did not significantly alter the dynamics of the flow, several trial runs must be conducted before finding the most appropriate configuration (number and position) of virtual atoms for a particular system.

Fig. 8 Computational domain of water confined between graphene sheets (dark blue spheres) and β-cristobalite (yellow spheres) wall with virtual particles (pink and light blue spheres). This figure has been adapted from ref. 133 with permission from American Chemical Society, copyright 2014.

4.2 Equilibrium molecular dynamics (EMD)

EMD simulations are used to calculate the friction coefficient at solid–liquid interfaces without relying on creating a non-equilibrium condition, either from flow-driven or pressure-driven setups. These methods do not resemble experimental setups as in NEMD, but are more reliable than NEMD simulations for systems with large values of slip, such as water flowing in carbon nanochannels,¹²⁷ where the velocity profiles are difficult to resolve. Notably, EMD methods often match NEMD calculations of slip in the low shear regime. The calculation of the friction factor using EMD is based on the fluctuation-dissipation theory, which leads to Green–Kubo-like expressions. The fluctuation-dissipation of particle forces, velocities, and their cross-correlations have been used to determine the friction factor employing four principal methods that are discussed herein.

Bocquet and Barrat¹¹⁸ proposed avoiding the arbitrariness of choosing the hydrodynamic boundary conditions by performing first-principles calculations. A “phenomenological” model of momentum transport based on the Navier–Stokes equations was formulated using a statistical mechanics approach, selecting L_S and the hydrodynamic distance to the wall (z_wall) as fitting parameters. Analytical expressions for the momentum density correlation function were obtained for perfectly flat and rough walls. MD simulations were performed to obtain the “exact” values of the momentum density correlation function. Equilibrium simulations of nano-confined atomic liquids interacting through purely repulsive potentials were conducted at constant temperature and varying the size of the channels. Through parameter fitting, Bocquet and Barrat¹¹⁸ calculated L_S and z_wall, and found that the analytical models matched the simulation results for repulsive walls with and without corrugation. However, confinements imposed by attractive walls were not correctly described by the phenomenological model due to the presence of slip-locking. Lastly, they conducted NEMD simulations of Couette flow to prove the effectiveness of the “phenomenological” model in predicting the velocity profile. The analytical model fitted through EMD simulations accurately matched the NEMD velocity profiles.

Given the success of their equilibrium calculations, Bocquet and Barrat¹¹⁸ formulated a method to calculate L_S and z_wall as equilibrium properties using a Green–Kubo-like approach. They employed linear response theory and the Mori–Zwanzig formalism separately to derive equilibrium coefficients based on the time-dependent correlation functions of the fluid. A perturbation Hamiltonian was chosen to generate a Poiseuille flow in the x-direction with a fictitious shear field [small gamma, Greek, dot above], where P_i,x is the momentum of particle i in the x-direction, and z₀ is the position at which the velocity profile vanishes while v_x(z) = [small gamma, Greek, dot above](z − z₀) is a first-order approximation of the tangential velocity. Applying linear response theory with a non-equilibrium friction force 〈F_x〉(t) as the response and H as the perturbation field, the following expression was developed:

where F_x represents the total force exerted on the solid wall by the liquid atoms in the x-direction; k_B, T, and A denote the Boltzmann constant, temperature, and the area of the wall, respectively; and the wall friction coefficient (λ) is:from which the slip length is obtained as L_S = η/λ.

According to Navier's friction model,⁴⁶F_x = −λAv_x(Δz), where F_x is the total tangential force exerted by the liquid on the wall and ν_x(Δz) is the tangential velocity at an equilibrium position Δz away from the wall; Bocquet and Barrat¹¹⁸ assumed that [small gamma, Greek, dot above] (z₀ − z_wall) = v_x(Δz) to find an expression for the friction coefficient. Petravic and Harrowell¹³⁴ indicated the incorrectness of such an assumption given that [small gamma, Greek, dot above] is an artificial constant field used to introduce a perturbation to the system and no physical correlation exists with F_x which is the friction force on the wall.

Petravic and Harrowell¹³⁴ addressed the equilibrium perturbation issue using Doll's equations of motion. They induced a disturbance within the system by generating a relative velocity, referred to as Δv_wall, between the confining walls, while also considering the constraints of a heterogeneous system in a boundary-driven flow context. The system's linear response was then evaluated in the context of a small Δv_wall, leading to the determination of a new friction coefficient.

where μ_i is determined when t → ∞ (statistically equilibrium stage). At this stage, μ₁(t) = μ₂(t) = μ (i = 1, 2 are the indices of two different confining walls). Eqn (9) bears similarity to eqn (7), with the distinction that Δv_wall is utilized in place of slip velocity to relate the friction force to velocity. Consequently, eqn (9) accounts for the entire thickness of the confined fluid and addresses the size-dependent friction coefficient observations noted in Petravic et al.¹³⁴ Finally, the L_S can be obtained as follows:where D represents the distance between the solid walls. eqn (10) can be simplified when considering identical walls.

Hansen et al.¹²⁵ emphasized that eqn (9) accounts for the friction of the whole system, including both the walls and the liquid, and highlighted the importance of separating the region affected by the wall from the bulk fluid to accurately determine the true wall friction. They expanded on Navier's foundational ideas, addressing the issue of wall friction by focusing on a thin layer of liquid close to the wall, see Fig. 9. The analysis considers the velocity profile of a liquid confined between two walls, separated by a distance L_y, where the liquid slab delimited at y = Δ is analyzed using Newton's second law:

where the mass of the liquid is m, while u_slab denotes the slab's center-of-mass velocity. The friction force resulting from the interactions between the wall and the slab is referred to as , whereas is the friction force due to the fluid-slab interactions. Lastly, F_e is an applied external force per unit mass.

Fig. 9 Sketch of the system used in the friction analysis in Hansen et al.¹²⁵ This figure has been redrawn from ref. 125 with permission from American Physical Society, copyright 2011.

The friction force was described by the following expression:

where ζ is a friction kernel, Δu = u_slab − u_wall, and represents a random force that has a zero mean and is uncorrelated with u_slab. In steady-state conditions, the friction forces are given as:where A = L_xL_z is the cross-section area of the system and ζ₀ is the zero-frequency friction coefficient.¹²⁵ Hansen et al.¹²⁵ calculated L_S using eqn (12–14), in combination with the Couette and Poiseuille flow solutions obtained with integral boundary conditions (finite liquid regions of width Δ). L_S calculated from such solutions matched the expected value of L_S = η/ξ₀ where ξ₀ = ζ₀/A in the limit when Δ → 0. The analytical solutions were compared with EMD simulations from which the friction coefficient was obtained from the Laplace transform of the velocity–velocity and force–velocity autocorrelation functions (ACFs) as:where the ACFs are defined as:

And

The Laplace transform of the friction kernel was obtained using a Maxwellian memory function for convenience as indicated in Hansen et al.,¹²⁵

Thus, the zero-frequency friction coefficient can be found by fitting B_i and κ_i in eqn (19) using data from EMD simulations as:

Hansen et al.¹²⁵ discovered that ζ₀ depended on Δ, requiring multiple trials to determine the appropriate Δ value. Thin slabs failed to capture the entire wall-slab interactions, while wider slabs may introduce unnecessary bulk particles. They noted that the friction coefficient was influenced by the channel width, particularly for channels with L_y ≤ 7σ, where σ denotes molecular diameter. A comparison was conducted between the L_S predicted using their EMD approach and NEMD calculations using Couette and Poiseuille flows, and a remarkable agreement between the two methods was found for flows with small shear rates.

Bocquet and Barrat¹³⁵ addressed the criticisms in the definition of their model reported in Bocquet et al.,¹¹⁸ suggesting that the sensitivity of their interfacial friction coefficient stemmed from the approach taken in handling system size and time limits extending to infinity, as derived from MD simulations. To reinforce the generality of their previous model, they introduced a new formulation for the Green–Kubo relationship for the friction coefficient λ, offering a more robust and fundamental approach grounded in the general Langevin equation. This formulation was applicable to both planar and cylindrically confined fluids. In their study, Bocquet and Barrat¹³⁵ defined a confined liquid system in which solid walls of large mass M are allowed to move in the tangential direction only. In the presence of solid–liquid friction, the fluctuations in the wall velocity U(t) are given by:

where the slip velocity is v_s(t) = U(t) − v_f(t) and v_f(t) is the fluid velocity, A is the wetted area, and δF(t) is a lateral fluctuating force. In the linear response regime, the slip velocity is related to the wall velocity as:where Ψ is a friction memory kernel related to the hydrodynamic shear modes in the fluid. After substituting eqn (21) into (20), the Langevin equation was Laplace-transformed and the force correlation function 〈F_w(t)F_w(0)〉 was found:

Eqn (23) is best handled in the Laplace space from which the friction coefficient was found when [small variant phi, Greek, tilde](s) is evaluated in the s → 0 limit, yielding:

which is the same expression previously reported by Bocquet et al.¹¹⁸ but this time using more general arguments and without the approximations involved in the first derivation.

Huang and Szlufarska¹²³ noted a significant concern in the discourse regarding friction coefficients derived from equilibrium calculations. They argued that the friction coefficient should be considered a local parameter rather than a bulk property. For instance, while solid–liquid friction exists in liquids moving through a carbon nanotube, achieving the thermodynamic limit in such a system, as proposed by Bocquet and Barrat,^118,135 is not feasible. Moreover, when dealing with heterogeneous surfaces or fluid mixtures in contact with a solid boundary, using a bulk property equation like eqn (24) fails to capture the localized variations at the interface where friction takes place. Bocquet et al.¹³⁵ applied the general Langevin equation along with a set of sum rules to the fluctuating velocity of a wall of large mass. In a new formulation, Huang and Szlufarska¹²³ applied a mechanical perturbation Hamiltonian to individual liquid particles at the solid–liquid interface, H = −xfe^iωt where x is the displacement of particles parallel to the solid walls and fe^iωt is an external drag force with frequency ω and time t. Being u_i the drift velocity of an interfacial particle moving parallel to the solid wall, and Fourier transforming the linear correlation function, the particle's drift velocity and mobility ϕ_i were obtained as:

Linear response theory was applied a second time using F_i as the force exerted by the wall on a single interfacial particle i yielding:

Now, by definition of the friction coefficient:

where the total friction coefficient can be obtained by summing all the contributing particles and normalizing it by the area of the interface:where the short-range nature of F_i allows to evaluate eqn (29) for any number of liquid particles, without compromising the interfacial aspects of the calculations. Huang and Szlufarska¹²³ observed difficulty in achieving a well-converged value of ϕ_i due to the particles not staying sufficient time near the wall and due to the sensitivity to the spatial definition of the interfacial region. This problem was solved by obtaining the right-hand side terms of eqn (29) through a Langevin formalism where single particles are analyzed. Using linear response theory and after many mathematical manipulations, they obtained:where the static friction factor [small eta, Greek, macron](0) can be used to obtain . Table 2 summarizes the EMD models for the calculation of L_S with key highlights.

Table 2 Overview of the theoretical models for computing friction factors via EMD

Source	Equations	Key highlights
Bocquet and Barrat^118		• Utilized both linear response theory and the Mori–Zwanzig formalism independently to derive equilibrium coefficients.
		• Introduced an artificial shear rate into the fluid using a Hamiltonian to simulate Poiseuille flow.
		• Analyzed the ACF of the friction factor focusing on the solid atoms.
		• Determined that the ACF should reach zero when t → ∞; and calculated λ at this point
Petravic and Harrowell^134		• Discovered that the ACF's integral does not reduce to zero over time but instead stabilizes at a constant value, implying a smooth decay of the force ACF.
		• Employed a method similar to Bocquet and Barrat^118 for calculating the friction factor, though they disagreed on the interpretation of the findings.
		• Identified that the friction coefficient varies with system size, indicating it is not an intrinsic interfacial property.
Hansen et al.^125		• Proposed isolating the region near the wall from the bulk to accurately determine wall friction.
		• Developed a dynamic analysis of a thin liquid slab adjacent to the wall, correlating the friction force with the slab velocity using a memory function.
		• Conducted the ACF analysis focusing on the interfacial liquid atoms.
Bocquet and Barrat^135		• Refined their earlier model to address criticisms regarding the generality of their Green–Kubo formulation.
		• Introduced a non-Markovian general Langevin framework to investigate perturbations in the wall velocity and slip behavior.
Huang and Szlufarska^123		• Applied linear response theory to a system of liquid particles subjected to perturbations by a Hamiltonian.
		• Assumed that particles interact independently, with interfacial interactions considered additive.
		• Conducted the ACF analysis on the liquid atoms.
		• The model indicated that several Langevin equations are needed to explore both wall velocity and slip behavior, indicating a linear relationship between these variables.

---

Bocquet and Barrat¹¹⁸ identified the deficiencies of integrating the force ACF in eqn (7) in the limit when the lag time goes to infinity. As an alternative, it was proposed to evaluate the integral only up to the point where the ACF reached its first zero. A year earlier, Español and Zúñiga¹³⁶ highlighted the issues with evaluating the friction factor integral from zero to arbitrary limits as indicated in Bocquet et al.¹¹⁸ The solid–liquid friction phenomenon was studied using Hamilton's equation with projection operators on a Brownian particle of infinite mass interacting with other particles. Complimentary EMD simulations were carried out by analyzing a fixed liquid particle (infinite mass particle) interacting with several other liquid particles to prove a correlation between the decay of the momentum ACF and the friction factor. The force ACF decreased fast and smoothly but the integral of such was rather noisy with a tendency to decay after long simulation times. No plateau of eqn (7) was found for long simulation times, as suggested when the thermodynamic limit was reached, but approximations of the friction factor could be extracted from shorter-time behaviors.

Español and Zúñiga¹³⁶ concluded that EMD calculations using Green–Kubo-like models are hindered by the order in which the thermodynamic limit (infinite number of particles) and the infinite time limit are taken since they do not commute. However, as the simulation systems get larger, it is expected that the friction coefficient calculations approach the obtained results in the limits discussed. Notably, several authors reported no issues with the evaluation of eqn (7).^{36,116,117,128,134} In these investigations, smooth time-dependent friction factors are reported with a plateau at which point the steady state friction factor is evaluated. Furthermore, consistency between EMD calculations using eqn (7) and NEMD has been reported.^116,117,128 Liang and Keblinski¹²⁸ obtained the friction factor of argon flowing between graphene surfaces and observed a steady friction factor plateau over a window of 12 ps analyzing data over 10 ns.

Tocci et al.³⁶ used ab initio MD and force field MD simulations to obtain the friction factor for graphene and hexagonal boron nitride in contact with water. By using both sources of atomic trajectories evaluated from 50 ps to 10 ns in a 1 ps time window, they obtained smooth friction coefficient integrals. Falk et al.¹¹⁷ evaluated the friction factor from force ACFs evaluated over 0.4 ns with a time window of 2 ps. It was indicated that at long timescales (typically nanoseconds), the integral vanishes due to the finite size of the system, but at intermediate times a plateau of the integral can be observed. Additionally, they did not observe confinement dependence on the friction factor. Contrariwise, Wei et al.¹¹⁶ used eqn (7) to determine the L_S in water confined between graphitic-carbon walls and observed a confinement effect on the friction coefficient after the viscosity was adjusted to the confinement level. The investigation by Harrowell and Petravic¹³⁴ focused on giving a better interpretation of eqn (7) and throughout their analysis, smooth time-dependent friction coefficients were observed. Contrarily, Huang and Szlufarska¹²³ observed rather fluctuating time-dependent friction coefficient graphs.

Arguments supporting and disproving the properties of the original friction coefficient expression derived by Bocquet and Barrat^118,135 can be found all over the literature. The vanishing behavior of the integral in eqn (7) has been consistently reported. Likewise, the confinement effect has been reported by some authors, but others did not capture that in their analysis. New analytical approaches and reinterpretations of the initial model have been proposed, but they have not been widely investigated. For example, only Kannam et al.^126,127 used the method, proposed by Hansen et al.,¹²⁵ to study the friction between liquids and graphite surfaces obtaining consistent results with NEMD simulations in the low shear rate limit. There is a notable debate surrounding the EMD analysis of hydrodynamics in nanoconfined liquids, and more comprehensive studies are necessary to reach definitive conclusions. The methodologies for analysis and simulation are not thoroughly detailed in existing literature, and the inconsistencies observed across various studies may stem from errors in postprocessing or data sampling during EMD simulations.

5 The hydrodynamics of nanoconfined flows

5.1 Hydrodynamic slip mechanisms and molecular origins

Thompson and Robbins^98,137 made significant contributions to the understanding of the stick-slip mechanisms of liquids moving past solid surfaces by examining this phenomenon from a thermodynamic standpoint rather than considering it as a hydrodynamic instability. Through NEMD simulations of Couette flow, they recorded the friction force, wall displacement, and structure factor over time as the wall velocity varied. Their findings revealed that increasing the binding strength between solid and liquid atoms led to solid-to-liquid transitions in the liquid particles at the interface. In some extreme cases, crystallization of interfacial liquid particles occurred, which was then disrupted by the high shear stresses present in the Couette flow. This interplay between solid–liquid binding and shear-induced disruption resulted in periodic phase transitions, thereby framing hydrodynamic slip within a thermodynamic context.

Lichter et al.¹³⁸ proposed that liquid molecules spend sufficient time near the wall to warrant a dynamic treatment of their molecular motion, based on the observed ordering of liquid particles in the direction perpendicular to the wall and the mass flux towards the solid. They developed a stochastic differential-difference equation for particles in the first adsorption layer near the wall, allowing for mass exchange between the bulk and interfacial particles. This approach was termed the variable-density Frenkel-Kontorova model (vdFK). The vdFK model was qualitatively successful in predicting the relationship between shear rate and L_S observed in NEMD simulations. Moreover, the model identified two distinct slip mechanisms observed in NEMD simulations: (1) slip caused by localized defect propagation, where particle exchange occurs between interfacial vacancies and the bulk, and (2) simultaneous slip of large liquid regions. At low shear rates, localized defects emerge within the liquid layer, with adjacent molecules quickly filling the resultant vacancies, as depicted in Fig. 10(a). This defect propagation is notably slow under low-shear conditions. In contrast, at high shear rates, the shear forces are sufficient to induce concurrent slip across large domains of the liquid layer, see Fig. 10(b). Fig. 10(c) illustrates the response of L_S to the applied force as modeled by the vdFK model, demonstrating that at low levels of force, L_S remains relatively constant. This stability is due to the sparse nature of molecular defects, which propagate slowly through the adsorbed liquid layer without significantly affecting the overall slip behavior. These defects do not cover a substantial area, thus minimally impacting the bulk liquid behavior. However, as the force increases, a sharp transition occurs due to the intensification of local defects. Ultimately, the system reaches a new plateau at higher applied forces, indicating that the liquid layer moves uniformly over the solid surface. Further increases in force do not significantly impact the L_S, suggesting a saturation of mobility mechanisms at the interface.

Fig. 10 Trajectories of liquid atoms predicted by the vdFK model, (a) at a low shear rate (defect propagation stage), and (b) at a high shear rate (translation of the entire first liquid layer or concurrent slip). The black solid dots represent the positions of liquid atoms at the interface. For an eye guide, the blue dots indicate the generating vacancies and filling the resultant vacancies. The vertical red solid lines illustrate the positions of seven solid atoms at the interface. (c) L_S as a function of the applied force calculated from the vdFK model. The L_Svs. force data were normalized to collapse the low-forcing data points onto a single curve for different ground states which denotes the ratio of liquid–liquid and solid–solid spacing. The arrows in panel (c) indicate the L_S for the corresponding molecular trajectories of the panel (a) and (b). This figure has been adapted from ref. 138 with permission from American Physical Society, copyright 2004.

Martini et al.¹³⁹ reported defect slip (like in the vdFK model) using low shear rate NEMD simulations. In this regime, liquid particles adjacent to the solid surface hop between equilibrium sites within a potential field generated by the solid following Arrhenius dynamics. At higher shear rates, they observed global slip, where the entire layer of liquid particles moves collectively. At the smallest wall velocity, atom movement was almost indiscernible, with minimal movement either upstream or downstream. As the wall velocity increased to intermediate levels (5 ms⁻¹ and 50 ms⁻¹), atoms displayed periods of stillness interspersed with sudden downstream shifts. At 50 ms⁻¹, the behavior began to show collective trends where groups of atoms might slip simultaneously, indicating the onset of a more coordinated movement. At the highest wall velocity simulated, a distinct global slip was observed where all atoms within the first adsorbed liquid layer move uniformly downstream, presenting parallel trajectories that indicate a cohesive and uniform motion over the solid wall. These observations align with the vdFK model. Additionally, Martini et al.¹³⁹ found a critical wall velocity for their system—a specific wall velocity that demarcates the transition from defect-driven slip to this observed global slip.

When molecular vacancies at the interface are widely spaced during liquid slip, the movement of a single atom from one equilibrium position to another happens independently. This independence allows for studying the dynamics of individual atoms and the application of transition-state theory.¹³⁹ Babu and Sathian¹²⁰ utilized Eyring's theory of reaction rates (transition-state theory), which models viscous flow as a chemical reaction where the primary process involves molecular hopping between equilibrium positions, see Fig. 11. In this context, liquid molecules must surpass an energy barrier created by neighboring molecules to reach a new equilibrium position. A comprehensive analytical model comprising six equations was developed, with shear viscosity and the friction coefficient being the primary outputs. NEMD simulations of water confined between graphene sheets and carbon nanotubes (CNTs) were conducted to directly compute the friction coefficient and estimate the activation energy of the liquid molecules—an essential input for the analytical model. The model's validity is contingent upon maintaining a low shear rate to ensure defect slip, as noted by Martini et al.¹³⁹ The friction coefficient predictions from the analytical model generally aligned with the numerical simulations, although there were instances of underestimation and overestimation for the various confinement levels studied. The authors indicated that they used different driving forces, in a Poiseuille flow configuration, and varied the channel dimensions. If this process is not performed carefully, markedly different shear rates would be produced if the driving force is kept constant while the channel size varies. Additionally, Babu and Sathian¹²⁰ found size-dependent friction coefficients for flow between graphene sheets but Falk et al.¹¹⁷ reported otherwise.

Fig. 11 Schematic of atom transition from one equilibrium position to another following Eyring's theory of reaction rates. The figure has been adapted from ref. 120 with permission from American Physical Society, copyright 2012.

Understanding friction forces at solid–liquid interfaces remains a significant challenge. For instance, Tocci et al.¹⁴⁰ investigated the friction of water at the interfaces of graphene and hexagonal boron nitride (h-BN) using ab initio molecular dynamics (AIMD) simulations. Notably, graphene and h-BN generated a similar interfacial water structure. Furthermore, the AIMD calculations revealed nearly identical contact angles for water droplets on graphene and h-BN sheets.¹⁴¹ Despite interfacial liquid structure and wettability similarities, the calculated friction coefficient on h-BN was approximately three times higher than on graphene. This significant difference was attributed to the greater corrugation of the energy landscape generated by h-BN, determined by the differences in the electronic structure of the two 2-D materials. To further investigate this phenomenon, Secchi et al.¹⁴² conducted experiments on water transport inside carbon nanotubes (CNTs) and boron nitride nanotubes (BNNTs). Their study revealed a significant radius-dependent slippage in CNTs, where water flowing through the nanotubes exhibits nearly frictionless interfaces, leading to exceptionally high flow rates. In contrast, BNNTs showed almost no slippage, despite their similar crystallography to CNTs. This difference highlights the influence of subtle atomic-scale interactions at the solid–liquid interface, suggesting a connection between hydrodynamic behavior and the electronic properties of the confining material.

Recently, Kavokine et al.¹⁴³ developed a quantum theory of the solid–liquid interface and introduced a new concept of quantum friction caused by the coupling of charge fluctuations in water to electronic excitations in the solid surface. In this theory, the authors argued that hydrodynamic friction arises not only from the static roughness of a solid surface (classical friction) but also from the interaction between water fluctuations and solid electronic excitations (quantum friction). Thus, this concept could be understood from the electronic contribution to the solid–liquid friction behavior. They also investigated water interactions with graphitic materials, where (i) graphene, a 2-D material, exhibited very low energy excitations at very small momenta (q ≤ 0. 02 Å⁻¹); suggesting that the electronic excitations in graphene are less likely to interact with water molecules over large distances, contributing negligible quantum friction compared to classical friction. Conversely, (ii) graphite exhibits unexpectedly high friction compared to graphene caused by the distinct electronic structure due to the coupling between its layers; this coupling leads to the emergence of low-energy plasmon modes in graphite, which are absent in single-layer graphene. These low-energy excitations, particularly the surface plasmon modes, strongly interact with water molecules at the interface. In graphite, the low-energy plasmon mode has a frequency of approximately 50 meV and is polarized perpendicularly to the layers. This mode has a flat dispersion over a range of momenta, meaning it can interact more effectively with the fluctuating electric fields of water molecules, particularly with the Debye mode of water. The strong interaction between the graphite plasmon modes and the water Debye mode leads to enhanced quantum friction at the interface. Thus, the overall friction at the graphite-water interface is higher than at the graphene-water interface. Furthermore, Bui et al.¹⁴⁴ applied a classical model that adjusts the dielectric properties of a solid using a simple model of charge density fluctuations in a carbon substrate. Their findings showed an increase in interfacial friction consistent with recent theories of quantum friction, with friction rising as the solid's dielectric spectrum features overlap with the librational and Debye modes of water.

5.2 Solid–liquid affinity characterized via wettability and liquid structuring effects on slip

In nanoconfined liquids, surface effects are predominant and one of the most significant is the solid–liquid affinity. A macroscopic outcome of such affinity can be characterized using the contact angle (surface wettability). From an experimental point of view, it is very difficult to change the wettability of a surface; however, MD simulations offer several options to do this, e.g., (i) modifying the solid–liquid atomic force field, (ii) manipulating the electrostatic interactions between solid and liquid particles, (iii) varying the surface atomic density, and (iv) modifying the simulation temperature. Although these simulations are limited to atomically smooth surfaces, important investigations have been conducted in this area.^145–147

Voronov et al.^96,97 used standard EMD simulations to determine the contact angle of a simple Lennard-Jones (LJ) fluid on a graphite-like solid (i.e., droplet wettability). A parametric analysis was conducted in which the solid–liquid LJ parameters were independently varied to assess their effects on the calculated contact angle. Increasing the value of the solid–liquid energy parameter (ε_sl) generated more hydrophilic surfaces and prompted liquid particles near the wall to mimic the wetted solid structure, as reported by Thompson et al.⁹⁸ A linear dependence of the contact angle on ε_sl was found. Alternatively, the LJ length parameter (σ_sl) produced changes in the surface energy landscape. Larger values of σ_sl mimicked smoother and more hydrophilic surfaces than smaller values of this parameter, similar observations are reported by Zhang et al.¹²⁴ Thus, two opposite trends were found depending on how the surface wettability is altered, and caution was advised for modeling slip surfaces. L_S increased as the contact angle increased when one modifies ε_sl; however, L_S decreased as the contact angle decreased when σ_sl was increased. Two different mechanisms are responsible for such behaviors, one is pertinent to a binding energy effect (ε_sl), and the other is relevant to the surface energy landscape granularity (σ_sl). Hydrophilic surfaces generated by smooth energy landscapes cause large slip; alternatively, hydrophilic surfaces generated by strong solid–liquid affinity led to small slip. Thus, liquids can slip over hydrophilic surfaces and hydrophobic surfaces can have minimum slip if the surface energy landscape allows for liquid particles to be trapped.

L _S is greatly influenced by the magnitude and type of the solid–liquid force field parameters, which define the interface affinity. Previous investigations correlated L_S at various interfaces to surface wettability. In a significant development, Huang et al.¹⁰³ made an important contribution by proposing a quasi-universal scaling relation that suggests that L_S is a function of wettability, L_S∼ (1 + cos θ)⁻², where θ denotes the contact angle, see Fig. 12(a). Ho et al.⁹⁹ challenged this relation by modifying the wettability of MgO through adjustments to its lattice constant, finding that L_S increased in more hydrophilic surfaces. In a related study, Wang et al.¹⁴⁸ used MD simulations and AFM experiments to determine the friction coefficient at various solid-water interfaces. Their findings revealed a significant limitation of using the contact angle alone to explain variations in friction coefficients at the nanoscale. Despite observing a similar contact angle, the friction coefficient increased 41 times as the surface charge increased from 0e to 0.36e. This rise in friction was attributed to localized potential energy fluctuations, which create additional energy barriers for water molecules, underscoring the limitation of wettability metrics to explain the friction coefficient and hydrodynamic slip in nanochannels. Wang et al.¹⁴⁹ further investigated the role of ordered water molecules at the solid–liquid interface of superhydrophilic surfaces using NEMD simulations. They observed that the formation of a hexagonal-like structure in the first water monolayer significantly reduced friction between the monolayer and bulk water above by decreasing the number of hydrogen bonds. The weakened hydrogen bonding led to smoother interlayer movement, thereby considerably reducing the overall friction at the interface. Supporting these observations, Xu et al.¹⁵⁰ conducted MD simulations in different polygonal carbon nanotubes (CNTs), demonstrating a similar frictional reduction due to the ordering of water molecules in the first monolayer.

Fig. 12 (a) Quasi-universal relationship, where L_S is a function of contact angle, proposed by Huang et al. The panel (a) has been adapted from ref. 103 with permission from American Physical Society, copyright 2008. L_S of Si and graphene-coated Si channels as a function of (b) the contact angle, and (c) the density depletion length δ. The (b) and (c) panels have been adapted from ref. 151 with permission from AIP publishing, copyright 2016.

These findings indicate that friction and hydrodynamic slip at the solid–liquid interface are affected by the liquid structuring at the interface rather than by wettability. To quantify this phenomenon, recent investigations introduced the concept of density depletion length (δ),^{93,103,104,151} which quantifies the presence (excess/deficit) of momentum-carrying liquid molecules at the interface (see Fig. 3c). The following equation has been reported to calculate δ:

where ρ_S and ρ_L represent the solid and liquid density distribution, respectively, with the superscript ‘b’ denoting a bulk value, which is characteristic of regions far from the interface. A lower δ value indicates a higher concentration and closer proximity of liquid particles to the solid surface, enhancing momentum transfer, while a higher δ value suggests fewer momentum carriers at the interface.

In addition to a quasi-universal relationship, Huang et al.¹⁰³ and Sendner et al.¹⁰⁴ reported that L_S correlates with the density depletion length as L_S∼ δ⁴. This scaling law is based on analyses using a mean-field theory model of wettability and a Green–Kubo-like model of slip, effectively explaining L_S behaviors across different models. Ramos-Alvarado et al.,¹⁵¹ in their series of EMD simulations on various Si nanochannels (bare Si (100) and Si (111)), and graphene-coated Si, noted that the quasi-universal relationship to θ only traced the data trends with limited fidelity and broke down for graphite-coated Si surfaces, as illustrated in Fig. 12(b). However, the scaling law of δ reliably quantified L_S across these diverse nanochannels, see Fig. 12(c).

In more detailed studies, Paniagua et al.⁹³ utilized a range of Lennard-Jones (LJ) parameters to model graphite-water interactions through EMD simulations. They reported that surface wettability was inadequate in characterizing L_S. Despite accurately controlling the surface wettability in their MD simulations (symbols inside the rectangular box in Fig. 13(a)), considerable variations in L_S were observed—26.87 nm, 42.61 nm, and 62.48 nm. These slip variations existed even when the contact angle, binding energy, and work of adhesion were similar across the three highlighted interfaces. This variability, driven by different friction coefficients, highlights the inadequacy of using wettability metrics alone to explain hydrodynamic slip behavior. Conversely, as depicted in Fig. 13(b), δ effectively captured the variations in L_S calculated via EMD across the different interface models, where wettability metrics fell short.

Fig. 13 L _S as a function of (a) the contact angle and (b) density depletion length. The (a) and (b) panels have been adapted from ref. 93 with permission from American Chemical Society, copyright 2020. (c) Shear-dependent L_S as a function of δ. The solid and dashed lines represent the averaged L_S at high and low shear rates. The (c) panel has been adapted from ref. 129 with permission from AIP publishing, copyright 2024. Each symbol represents a different interface model.

Furthermore, Paniagua et al.⁹³ noted that while the wettability scaling law could generally describe the behavior across most interface conditions, it failed in extreme hydrophobic or hydrophilic scenarios. To address this, they proposed an empirical exponential function (L_S∼ e^Bδ) that could effectively trace L_S across all graphite-water interface models. Corroborating these findings, Shuvo et al.¹²⁹ conducted NEMD simulations in shear-driven flows within graphite nanochannels, confirming that the exponential function of δ also accurately describes the behavior of L_S under different shear conditions, as illustrated in Fig. 13(c).

5.3 Shear rate effect on hydrodynamic slip

Shear rate is crucial in defining the boundary condition in nanoconfined liquids as the shear force on liquid particles competes with the solid–liquid binding and liquid–liquid cohesive forces. Thompson and Troian¹⁰⁰ carried out NEMD simulations of an LJ liquid, varying parameters like ε_sl, σ_sl, and solid density using different shear rates in a Couette flow model. They observed that L_S remained constant over a certain shear rate range but exhibited rapid growth beyond a critical value γ_c. Despite the variation in shear rates, the shear viscosity showed no significant change, indicating Newtonian behavior. A universal boundary condition was proposed: , where L⁰_S represents the low shear rate L_S limit and α is a fitting parameter, suggesting that the Navier slip condition is just a specific case of a broader relationship.

The rapid increase in L_S at a critical shear rate was similarly observed in studies by Voronov et al.^96,97 and Chen et al.³³ in shear-driven MD simulations. Kannam et al.¹²⁶ investigated hydrodynamic slip for both Poiseuille and Couette flows of graphite–argon, and graphite–methene systems, reporting an exponential (unbounded) growth of the L_S in both flow types. The authors did not address the seemingly infinite growth of the L_S when factors such as the wall friction coefficient and the fluid's viscosity pose a physical limit to solid–liquid friction. Wagemann et al.¹⁵² expanded on this by examining L_S of water within graphene nanochannels, particularly focusing on the crystallographic directions—zig-zag and arm-chair. Their observations indicated an unbounded growth of L_S in both directions, see Fig. 14(a). The authors calculated the wall friction coefficient as a function of shear rate and found that at low shear rates, the friction coefficient remained constant, indicating a stable interaction between the fluid and the solid surface. However, at high shear rates, a rapid reduction in the friction coefficient was observed, suggesting an unbounded growth of L_S. Notably, the authors did not investigate the rheological properties of the liquid, which are crucial because L_S is a function of both fluid viscosity and wall friction coefficient. Recently, Li et al.¹⁵³ investigated the rheology of water in nanoconfined graphite walls and suggested a shear thinning effect at high shear rates.

Fig. 14 Shear effect on slip reported through the years: (a) unbounded growth of the L_S, (this figure has been adapted from ref. 152 with permission from Royal Society of Chemistry, copyright 2017); (b) bounded growth of the L_S, (this figure has been adapted from ref. 106 with permission from American Physical Society, copyright 2008); and (c) reduction of the L_S, (this figure has been adapted from ref. 154 with permission from American Physical Society, copyright 2011).

Conversely, Martini et al.¹³⁹ observed different regimes of slip featuring a bounded growth of L_S after a given critical shear rate, see Fig. 14(b), corroborating their molecular mechanism of slip theory. In the literature concerning the MD modeling of droplet wettability, it was reported that keeping the solid atoms rigid not only allowed to significantly reduce the computing times but also a negligible variation of the contact angle was observed compared with flexible wall models.³³ Thus, several early contributions took a similar approach for their NEMD simulations of slip.^96,97 Martini et al.¹⁰⁶ hypothesized that the unbounded growth of L_S with increasing shear rate, observed in previous simulations, was due to using rigid wall atoms, which overlooked momentum transfer between solid and liquid particles. To validate this, they conducted NEMD simulations of Couette flow with both fixed and flexible wall atoms. The results confirmed that rigid walls lead to unbounded L_S growth at high shear rates, whereas flexible walls exhibited a constant L_S beyond a certain shear rate threshold.

Pahlavan and Freund¹⁵⁴ suggested reevaluating the high shear rate limit in NEMD simulations by decoupling the effects of the wall and thermostating approaches. Their findings indicated that the solid–liquid vibrational frequency mismatch had a negligible effect on L_S, while a reduction in L_S was attributed to the local temperature rise caused by an increasing number of solid–liquid collisions at high shear rates, see Fig. 14(c). Furthermore, the reduction of L_S was also reported by Ramos-Alvarado et al.¹⁵⁵ through NEMD simulation of both Couette and Poiseuille flow.

Further detailed investigations by Shuvo et al.¹²⁹ using different interface models of graphite-water interfaces under different shear conditions showed a bounded growth of the L_S at high shear rates, see Fig. 15(a), aligning with Martini's findings using the flexible wall model. They explored the rheology of water and the wall friction coefficient to understand the bimodal response of L_S under varying shear conditions. They discovered that both viscosity and friction coefficient decreased at high shear rates but remained constant at lower shear rates, see Fig. 15(b). During the transition from low to high shear rates (LSR to HSR), the friction coefficient decreased more rapidly than the shear viscosity, until reaching a new equilibrium. As a result, the L_S was higher and constant at higher shear rates.

Fig. 15 (a) The bimodal response of L_S under different shear conditions. (b) Normalized viscosity and friction coefficient in shear-driven flow. The figure has been adapted from ref. 129 with permission from AIP publishing, copyright 2024. (a) was redrawn with sigmoid fits using the data reported by Shuvo et al.¹²⁹

As discussed in section 4.1, NEMD simulations are pivotal in calculating the transport properties of molecular systems by mimicking experimental setups. However, achieving a good signal-to-noise ratio necessitates applying an external perturbation that is significantly larger than those typically used in experiments. This approach helps overcome the limitations of short simulation timescales and smaller length scales compared to those in experimental settings, but it also presents a challenge for directly validating the simulation results. Addressing this limitation, Maffioli et al.¹⁵⁶ developed the TTCF4LAMMPS technique, which combines direct NEMD simulations with the Transient Time Correlation Function (TTCF). This integration facilitates the exploration of fluid responses at shear rates achievable in experiments. TTCF relies on the correlation between the initial rate of energy dissipation and the response of any phase variable following an external perturbation, described mathematically as:

where B(t) represents an arbitrary dynamic variable of interest, and Ω denotes the dissipation function related to the system energy changes due to external perturbations. 〈Ω(0)B(s)〉 represents the cross-correlation between the initial dissipation and the variable at time s. In the TTCF methodology, mother and daughter trajectories are essential for analyzing system properties. Initially, a mother trajectory is established through EMD simulations, allowing the system to evolve under equilibrium conditions to provide a statistical baseline. From this, several daughter trajectories are generated at varied intervals, each inheriting initial states from the mother trajectory but experiencing specific external perturbations. These daughter trajectories are crucial for examining the system's response to these perturbations, with their transient responses averaged to determine the desired transport properties of the system with a good signal-to-noise ratio at low (realistic) shear rates.

Despite the advantages of TTCF over traditional NEMD in terms of accessing experimentally relevant shear rates, its adoption remains limited, possibly due to the complexity of its implementation and the high computational demands associated with it.

6 Summary and outlook

Research into fluid dynamics at the nanoscale has revealed notable deviations from continuum fluid behavior, particularly regarding the phenomenologically reported no-slip condition at the solid–liquid interface. At such scales, surface effects—such as roughness, wettability, and molecular interactions—become increasingly significant due to the dimensions being on the order of molecular mean free paths. This leads to unique properties like altered viscosity and density profiles near the interface, complicating the understanding of flow dynamics.

In Section 3, experimental methodologies including SFA, AFM, and μ-PIV were introduced for the direct measurement of the hydrodynamic slip length. SFA and AFM quantify the viscous drag force, which is balanced by the restorative force of the AFM cantilever based on the separation between the surface and the AFM tip. μ-PIV profiles the velocity of liquid near a solid wall by tracking particle movement. Additionally, new approaches using microchannel resonators and hybrid nanochannels with varying flow resistances were presented. These techniques are highly valuable for studying nanoscale hydrodynamic phenomena, but they also have limitations in providing detailed information on the interfacial interactions necessary to fully understand the origins of hydrodynamic slip. In Section 3.5, SFG and XRR were introduced as typical interfacial analysis techniques to study water molecules at solid surfaces. However, it remains unclear how the structural and chemical interactions of interfacial water molecules influence hydrodynamic slip behavior. These challenges have steered the field towards computational methods, where NEMD and EMD simulations play pivotal roles. NEMD simulations, which simulate shear rate effects in a manner akin to experimental setups, often require large velocity gradients to mitigate statistical noise, which can lead to unphysical conditions. Alternatively, EMD simulations focus on tracking the linear response of systems, providing a more reliable means of computing transport properties, such as interfacial friction coefficients. However, the application of Green–Kubo relations in these simulations must be handled with care to avoid introducing non-physical parameters that could skew the interpretation of nanoconfined flow characteristics.

The effect of increasing shear rate on hydrodynamic slip varies widely depending on simulation conditions, with reports of unbounded, bounded, and even reduced L_S. The unbounded L_S is often an artifact of simulations that neglect momentum transfer between solid and liquid particles, while reduced L_S is linked to local temperature increases due to frequent solid–liquid collisions at high shear rates. The bounded growth of L_S has been supported by MD simulations with flexible solid wall models where momentum transfer between solid and liquid is allowed, experiments, and theoretical models like the vdFK. Furthermore, EMD-computed L_S matches with NEMD results at low shear rates. A significant limitation of traditional NEMD simulations is their reliance on very large velocity gradients. The recent development of the TCFF4LAMMPS technique enables the generation of numerical data with high signal-to-noise ratios at velocity gradients that are accessible in experiments, thereby facilitating the validation of computational models with experimental findings.

Significant advancements have been made in MD simulations unveiling atomistic details of hydrodynamic slip behavior at the solid–liquid interface, yet several fundamental questions remain unresolved. The challenge lies in developing a physics-informed boundary condition that accounts for complex interfacial interactions, including surface chemistry and interfacial liquid structuring. Current research indicates that a liquid structuring parameter may play more significant roles in hydrodynamic slip than wettability metrics, a hypothesis requiring verification. Inter-particle interactions in the liquid near the solid surface cannot be fully modeled with the parameters determined from the bulk phase properties. The solid–liquid interaction right at the interface can influence the intermolecular interaction in the next layer, which will propagate further into the liquid phase. How fast or slow this interaction decays with the distance from the surface could be another important parameter that governs (at least affects) how effectively the momentum will be transferred from the bulk liquid to the solid surface. Although indirect, such interactions could be extracted from advanced characterization methods that are sensitive to structural order or density change in the liquid phase in proximity to the solid surface. By integrating the solid–liquid interaction parameters with the boundary conditions, MD simulations will be able to predict and explain hydrodynamic slip behavior at the liquid/solid interface.

Author contributions

Abdul Aziz Shuvo: data curation (experimental and numerical part, equal), formal analysis (experimental and numerical part, equal), investigation (experimental and numerical part, equal), writing – original draft (experimental and numerical part, equal). Luis E. Paniagua-Guerra: data curation (experimental and numerical part, equal), formal analysis (experimental and numerical part, equal), investigation (experimental and numerical part, equal), writing – original draft (experimental and numerical part, equal), writing – review & editing (equal). Juseok Choi: data curation (experimental part, equal), formal analysis (experimental part, equal), investigation (experimental part, equal), writing – original draft (experimental part, equal). Seong H. Kim: data curation (experimental part, equal), formal analysis (experimental part, equal), investigation (experimental part, equal), writing – review & editing (equal), funding acquisition (equal), supervision (equal). Bladimir Ramos-Alvarado: data curation (experimental and numerical part, equal), formal analysis (experimental and numerical part, equal), investigation (experimental and numerical part, equal), writing – original draft (experimental and numerical part, equal), writing – review & editing (equal), funding acquisition (equal), supervision (equal).

Data availability

No primary research results, software or code have been included and no new data were generated or analyzed as part of this review.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Science Foundation, USA (Award number: 2241730), and by the Materials Research Institute at Penn State University through their Seed Grant Program.

References

1. C. L. Rice and R. Whitehead, J. Phys. Chem., 1965, 69, 4017–4024.
2. J. C. T. Eijkel and A. van den Berg, Microfluid. Nanofluid., 2005, 1, 249–267.
3. L. Bocquet and E. Charlaix, Chem. Soc. Rev., 2010, 39, 1073–1095.
4. L. J. Guo, X. Cheng and C.-F. Chou, Nano Lett., 2004, 4, 69–73.
5. W. Sun, P. Qin, H. Gao, G. Li and K. Jiao, Biosens. Bioelectron., 2010, 25, 1264–1270.
6. V. Sunkara, B. J. Hong and J. W. Park, Biosens. Bioelectron., 2007, 22, 1532–1537.
7. P. M. Sinha, G. Valco, S. Sharma, X. Liu and M. Ferrari, Nanotechnology, 2004, 15, S585.
8. R. Lopez-Salazar, S. Camacho-Leon, L. Olivares-Quiroz and J. Hernandez, Procedia Technol., 2012, 3, 334–341.
9. W.-H. Lee, C.-Y. Loo, D. Traini and P. M. Young, Expert Opin. Drug Delivery, 2015, 12, 1009–1026.
10. Z. Mazibuko, Y. E. Choonara, P. Kumar, L. C. Du Toit, G. Modi, D. Naidoo and V. Pillay, J. Pharm. Sci., 2015, 104, 1213–1229.
11. G. Wang, W. Mao, R. Byler, K. Patel, C. Henegar, A. Alexeev and T. Sulchek, PLoS One, 2013, 8, e75901.
12. V. Soum, S. Park, A. I. Brilian, O.-S. Kwon and K. Shin, Micromachines, 2019, 10, 516.
13. P. Yager, T. Edwards, E. Fu, K. Helton, K. Nelson, M. R. Tam and B. H. Weigl, Nature, 2006, 442, 412–418.
14. M. L. Kovarik and S. C. Jacobson, Anal. Chem., 2009, 81, 7133–7140.
15. H. Daiguji, P. Yang, A. J. Szeri and A. Majumdar, Nano Lett., 2004, 4, 2315–2321.
16. A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell and L. Bocquet, Nature, 2013, 494, 455–458.
17. B. E. Logan and M. Elimelech, Nature, 2012, 488, 313–319.
18. K.-H. Paik, Y. Liu, V. Tabard-Cossa, M. J. Waugh, D. E. Huber, J. Provine, R. T. Howe, R. W. Dutton and R. W. Davis, ACS Nano, 2012, 6, 6767–6775.
19. R. Fan, M. Yue, R. Karnik, A. Majumdar and P. Yang, Phys. Rev. Lett., 2005, 95, 086607.
20. R. Karnik, R. Fan, M. Yue, D. Li, P. Yang and A. Majumdar, Nano Lett., 2005, 5, 943–948.
21. D. Constantin and Z. S. Siwy, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 041202.
22. I. Vlassiouk and Z. S. Siwy, Nano Lett., 2007, 7, 552–556.
23. R. Karnik, C. Duan, K. Castelino, H. Daiguji and A. Majumdar, Nano Lett., 2007, 7, 547–551.
24. G. Nguyen and Z. Siwy, Biophys. J., 2010, 98, 602a.
25. R. R. Nair, H. A. Wu, P. N. Jayaram, I. V. Grigorieva and A. K. Geim, Science, 2012, 335, 442–444.
26. D. Mijatovic, J. C. T. Eijkel and A. van den Berg, Lab Chip, 2005, 5, 492.
27. B. Radha, A. Esfandiar, F. C. Wang, A. P. Rooney, K. Gopinadhan, A. Keerthi, A. Mishchenko, A. Janardanan, P. Blake, L. Fumagalli, M. Lozada-Hidalgo, S. Garaj, S. J. Haigh, I. V. Grigorieva, H. A. Wu and A. K. Geim, Nature, 2016, 538, 222–225.
28. D. Cohen-Tanugi and J. C. Grossman, Nano Lett., 2012, 12, 3602–3608.
29. D. Konatham, J. Yu, T. A. Ho and A. Striolo, Langmuir, 2013, 29, 11884–11897.
30. A. Aghigh, V. Alizadeh, H. Y. Wong, Md. S. Islam, N. Amin and M. Zaman, Desalination, 2015, 365, 389–397.
31. K. A. Mahmoud, B. Mansoor, A. Mansour and M. Khraisheh, Desalination, 2015, 356, 208–225.
32. S. Rikhtehgaran and A. Lohrasebi, Desalination, 2015, 365, 176–181.
33. Y. Chen, D. Li, K. Jiang, J. Yang, X. Wang and Y. Wang, J. Chem. Phys., 2006, 125, 084702.
34. M. Sega, M. Sbragaglia, L. Biferale and S. Succi, Soft Matter, 2013, 9, 8526.
35. W. Song, Y. Tang, C. Qian, B. J. Kim, Y. Liao and D.-G. Yu, Innovation, 2023, 4, 100381.
36. G. Tocci, L. Joly and A. Michaelides, Nano Lett., 2014, 14, 6872–6877.
37. J. K. Holt, H. G. Park, Y. Wang, M. Stadermann, A. B. Artyukhin, C. P. Grigoropoulos, A. Noy and O. Bakajin, Science, 2006, 312, 1034–1037.
38. M. Majumder, N. Chopra, R. Andrews and B. J. Hinds, Nature, 2005, 438, 44–44.
39. G. Hummer, J. C. Rasaiah and J. P. Noworyta, Nature, 2001, 414, 188–190.
40. J. A. Thomas and A. J. H. McGaughey, Phys. Rev. Lett., 2009, 102, 184502.
41. K. Wu, Z. Chen, J. Li, X. Li, J. Xu and X. Dong, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 3358–3363.
42. J.-J. Shu, J. B. M. Teo and W. K. Chan, Soft Matter, 2016, 12, 8388–8397.
43. E. Secchi, S. Marbach, A. Niguès, D. Stein, A. Siria and L. Bocquet, Nature, 2016, 537, 210–213.
44. M. Ma, F. Grey, L. Shen, M. Urbakh, S. Wu, J. Z. Liu, Y. Liu and Q. Zheng, Nat. Nanotechnol., 2015, 10, 692–695.
45. D. C. Tretheway and C. D. Meinhart, Phys. Fluids, 2002, 14, L9–L12.
46. P. M. Navier, Mem. Acad. Sci., 1823, 6, 389–440.
47. J. N. Israelachvili and G. E. Adams, Nature, 1976, 262, 774–776.
48. O. I. Vinogradova, J. Colloid Interface Sci., 1995, 169, 306–312.
49. D. Y. C. Chan and R. G. Horn, J. Chem. Phys., 1985, 83, 5311–5324.
50. G. Luengo, F.-J. Schmitt, R. Hill and J. Israelachvili, Macromolecules, 1997, 30, 2482–2494.
51. J. Baudry, E. Charlaix, A. Tonck and D. Mazuyer, Langmuir, 2001, 17, 5232–5236.
52. Y. Zhu and S. Granick, Phys. Rev. Lett., 2001, 87, 096104.
53. C. Cottin-Bizonne, B. Cross, A. Steinberger and E. Charlaix, Phys. Rev. Lett., 2005, 94, 056102.
54. C. Cottin-Bizonne, A. Steinberger, B. Cross, O. Raccurt and E. Charlaix, Langmuir, 2008, 24, 1165–1172.
55. V. S. J. Craig, C. Neto and D. R. M. Williams, Phys. Rev. Lett., 2001, 87, 054504.
56. O. I. Vinogradova, H.-J. Butt, G. E. Yakubov and F. Feuillebois, Rev. Sci. Instrum., 2001, 72, 2330–2339.
57. O. I. Vinogradova and G. E. Yakubov, Langmuir, 2003, 19, 1227–1234.
58. E. Bonaccurso, M. Kappl and H.-J. Butt, Phys. Rev. Lett., 2002, 88, 076103.
59. C. D. F. Honig and W. A. Ducker, Phys. Rev. Lett., 2007, 98, 028305.
60. B. Bhushan, Y. Wang and A. Maali, Langmuir, 2009, 25, 8117–8121.
61. A. Maali, C. Hurth, T. Cohen-Bouhacina, G. Couturier and J.-P. Aimé, Appl. Phys. Lett., 2006, 88, 163504.
62. A. Maali, T. Cohen-Bouhacina and H. Kellay, Appl. Phys. Lett., 2008, 92, 053101.
63. H. Li, Z. Xu, C. Ma and M. Ma, Nanoscale, 2022, 14, 14636–14644.
64. A. Maali, Y. Wang and B. Bhushan, Langmuir, 2009, 25, 12002–12005.
65. C. Zhang, X. Wang, J. Jin, L. Li and J. D. Miller, Colloids Interfaces, 2021, 5, 44.
66. C. D. F. Honig and W. A. Ducker, J. Phys. Chem. C, 2007, 111, 16300–16312.
67. O. I. Vinogradova, Langmuir, 1995, 11, 2213–2220.
68. H. Ishida, H. Teshima, Q.-Y. Li and K. Takahashi, Int. J. Thermofluids, 2024, 22, 100634.
69. D. Lumma, A. Best, A. Gansen, F. Feuillebois, J. O. Rädler and O. I. Vinogradova, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 056313.
70. J. Ou and J. P. Rothstein, Phys. Fluids, 2005, 17, 103606.
71. P. Joseph and P. Tabeling, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 035303.
72. J. F. Collis, S. Olcum, D. Chakraborty, S. R. Manalis and J. E. Sader, Nano Lett., 2021, 21, 4959–4965.
73. Q. Xie, M. A. Alibakhshi, S. Jiao, Z. Xu, M. Hempel, J. Kong, H. G. Park and C. Duan, Nat. Nanotechnol., 2018, 13, 238–245.
74. B. Radha, A. Esfandiar, F. C. Wang, A. P. Rooney, K. Gopinadhan, A. Keerthi, A. Mishchenko, A. Janardanan, P. Blake, L. Fumagalli, M. Lozada-Hidalgo, S. Garaj, S. J. Haigh, I. V. Grigorieva, H. A. Wu and A. K. Geim, Nature, 2016, 538, 222–225.
75. J. Ping and A. T. C. Johnson, Appl. Phys. Lett., 2016, 109, 013103.
76. A. Ambrosi, C. K. Chua, A. Bonanni and M. Pumera, Chem. Rev., 2014, 114, 7150–7188.
77. V. P. Zhdanov and B. Kasemo, Chem. Phys. Lett., 2011, 513, 124–126.
78. Q. Du, E. Freysz and Y. R. Shen, Phys. Rev. Lett., 1994, 72, 238–241.
79. Y. R. Shen and V. Ostroverkhov, Chem. Rev., 2006, 106, 1140–1154.
80. F. Wei, S. Urashima, S. Nihonyanagi and T. Tahara, J. Am. Chem. Soc., 2023, 145, 8833–8846.
81. A. T. Celebi, M. Barisik and A. Beskok, J. Chem. Phys., 2017, 147, 164311.
82. A. T. Celebi, M. Barisik and A. Beskok, Microfluid. Nanofluid., 2018, 22, 7.
83. X. Geng, M. Yu, W. Zhang, Q. Liu, X. Yu and Y. Lu, Sci. Rep., 2019, 9, 18957.
84. Y. Xie, L. Fu, T. Niehaus and L. Joly, Phys. Rev. Lett., 2020, 125, 014501.
85. M. Rezaei, A. R. Azimian and A. R. Pishevar, Phys. Chem. Chem. Phys., 2018, 20, 30365–30375.
86. Y. Wang, F. Tang, X. Yu, T. Ohto, Y. Nagata and M. Bonn, Angew. Chem., Int. Ed., 2024, 63, e202319503.
87. M. Masuduzzaman, C. Bakli, M. Barisik and B. Kim, Small, 2024, 2404397.
88. M. Masuduzzaman and B. Kim, Phys. Fluids, 2024, 36, 062003.
89. M. Masuduzzaman and B. H. Kim, Langmuir, 2022, 38, 7244–7255.
90. P. Ma, Y. Liu, X. Sang, J. Tan, S. Ye, L. Ma and Y. Tian, J. Colloid Interface Sci., 2022, 626, 324–333.
91. M. Mezger, H. Reichert, S. Schöder, J. Okasinski, H. Schröder, H. Dosch, D. Palms, J. Ralston and V. Honkimäki, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 18401–18404.
92. L. Wang, C. Zhao, M. H. G. Duits, F. Mugele and I. Siretanu, Sens. Actuators, B, 2015, 210, 649–655.
93. L. E. Paniagua-Guerra, C. U. Gonzalez-Valle and B. Ramos-Alvarado, Langmuir, 2020, 36, 14772–14781.
94. M. Mezger, B. M. Ocko, H. Reichert and M. Deutsch, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 3733–3737.
95. K. E. Karim, M. Barisik, C. Bakli and B. H. Kim, Phys. Chem. Chem. Phys., 2024, 26, 19069–19082.
96. R. S. Voronov, D. V. Papavassiliou and L. L. Lee, J. Chem. Phys., 2006, 124, 204701.
97. R. S. Voronov, D. V. Papavassiliou and L. L. Lee, Chem. Phys. Lett., 2007, 441, 273–276.
98. P. A. Thompson and M. O. Robbins, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 41, 6830–6837.
99. T. A. Ho, D. V. Papavassiliou, L. L. Lee and A. Striolo, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 16170–16175.
100. P. A. Thompson and S. M. Troian, Nature, 1997, 389, 360–362.
101. H. Zhang, Z. Zhang, Y. Zheng and H. Ye, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 066303.
102. R. Khare, P. Keblinski and A. Yethiraj, Int. J. Heat Mass Transfer, 2006, 49, 3401–3407.
103. D. M. Huang, C. Sendner, D. Horinek, R. R. Netz and L. Bocquet, Phys. Rev. Lett., 2008, 101, 226101.
104. C. Sendner, D. Horinek, L. Bocquet and R. R. Netz, Langmuir, 2009, 25, 10768–10781.
105. J. Xu and Y. Li, Int. J. Heat Mass Transfer, 2007, 50, 2571–2581.
106. A. Martini, H.-Y. Hsu, N. A. Patankar and S. Lichter, Phys. Rev. Lett., 2008, 100, 206001.
107. C. Liu and Z. Li, AIP Adv., 2011, 1, 032108.
108. C. Liu and Z. Li, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 036302.
109. C. Liu and Z. Li, J. Chem. Phys., 2010, 132, 024507.
110. G. Nagayama and P. Cheng, Int. J. Heat Mass Transfer, 2004, 47, 501–513.
111. G. Nagayama, T. Tsuruta and P. Cheng, Int. J. Heat Mass Transfer, 2006, 49, 4437–4443.
112. S. Ge, Y. Gu and M. Chen, Mol. Phys., 2015, 113, 703–710.
113. Z. Zhang, H. Zhang and H. Ye, Appl. Phys. Lett., 2009, 95, 154101.
114. Z. Li, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 026312.
115. A. E. Giannakopoulos, F. Sofos, T. E. Karakasidis and A. Liakopoulos, Int. J. Heat Mass Transfer, 2012, 55, 5087–5092.
116. N. Wei, X. Peng and Z. Xu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 012113.
117. K. Falk, F. Sedlmeier, L. Joly, R. R. Netz and L. Bocquet, Nano Lett., 2010, 10, 4067–4073.
118. L. Bocquet and J.-L. Barrat, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1994, 49, 3079–3092.
119. J. A. Thomas and A. J. H. McGaughey, Nano Lett., 2008, 8, 2788–2793.
120. J. S. Babu and S. P. Sathian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051205.
121. J.-L. Barrat and L. Bocquet, Faraday Discuss., 1999, 112, 119–128.
122. J.-L. Barrat and L. Bocquet, Phys. Rev. Lett., 1999, 82, 4671–4674.
123. K. Huang and I. Szlufarska, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 032119.
124. H. Zhang, Z. Zhang and H. Ye, Microfluid. Nanofluid., 2012, 12, 107–115.
125. J. S. Hansen, B. D. Todd and P. J. Daivis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 016313.
126. S. K. Kannam, B. D. Todd, J. S. Hansen and P. J. Daivis, J. Chem. Phys., 2011, 135, 144701.
127. S. K. Kannam, B. D. Todd, J. S. Hansen and P. J. Daivis, J. Chem. Phys., 2012, 136, 024705.
128. Z. Liang and P. Keblinski, J. Chem. Phys., 2015, 142, 134701.
129. A. A. Shuvo, L. E. Paniagua-Guerra, X. Yang and B. Ramos-Alvarado, J. Chem. Phys., 2024, 160, 194704.
130. S. Bernardi, B. D. Todd and D. J. Searles, J. Chem. Phys., 2010, 132, 244706.
131. X. Yong and L. T. Zhang, J. Chem. Phys., 2013, 138, 084503.
132. S. De Luca, B. D. Todd, J. S. Hansen and P. J. Daivis, J. Chem. Phys., 2014, 140, 054502.
133. S. De Luca, B. D. Todd, J. S. Hansen and P. J. Daivis, Langmuir, 2014, 30, 3095–3109.
134. J. Petravic and P. Harrowell, J. Chem. Phys., 2007, 127, 174706.
135. L. Bocquet and J.-L. Barrat, J. Chem. Phys., 2013, 139, 044704.
136. P. Español and I. Zúñiga, J. Chem. Phys., 1993, 98, 574–580.
137. P. A. Thompson and M. O. Robbins, Science, 1990, 250, 792–794.
138. S. Lichter, A. Roxin and S. Mandre, Phys. Rev. Lett., 2004, 93, 086001.
139. A. Martini, A. Roxin, R. Q. Snurr, Q. Wang and S. Lichter, J. Fluid Mech., 2008, 600, 257–269.
140. G. Tocci, L. Joly and A. Michaelides, Nano Lett., 2014, 14, 6872–6877.
141. H. Li and X. C. Zeng, ACS Nano, 2012, 6, 2401–2409.
142. E. Secchi, S. Marbach, A. Niguès, D. Stein, A. Siria and L. Bocquet, Nature, 2016, 537, 210–213.
143. N. Kavokine, M. L. Bocquet and L. Bocquet, Nature, 2022, 602, 84–90.
144. A. T. Bui, F. L. Thiemann, A. Michaelides and S. J. Cox, Nano Lett., 2023, 23, 580–587.
145. F. Taherian, V. Marcon, N. F. A. van der Vegt and F. Leroy, Langmuir, 2013, 29, 1457–1465.
146. C.-J. Shih, Q. H. Wang, S. Lin, K.-C. Park, Z. Jin, M. S. Strano and D. Blankschtein, Phys. Rev. Lett., 2012, 109, 176101.
147. T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu and P. Koumoutsakos, J. Phys. Chem. B, 2003, 107, 1345–1352.
148. C. Wang, H. Yang, X. Wang, C. Qi, M. Qu, N. Sheng, R. Wan, Y. Tu and G. Shi, Commun. Chem., 2020, 3, 27.
149. C. Wang, B. Wen, Y. Tu, R. Wan and H. Fang, J. Phys. Chem. C, 2015, 119, 11679–11684.
150. X. Xu, Z. Li, Y. Zhang, C. Wang, J. Zhao and N. Wei, Carbon, 2024, 228, 119402.
151. B. Ramos-Alvarado, S. Kumar and G. P. Peterson, Appl. Phys. Lett., 2016, 108, 074105.
152. E. Wagemann, E. Oyarzua, J. H. Walther and H. A. Zambrano, Phys. Chem. Chem. Phys., 2017, 19, 8646–8652.
153. F. Li, I. A. Korotkin and S. A. Karabasov, Langmuir, 2020, 36, 5633–5646.
154. A. Pahlavan and J. B. Freund, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 021602.
155. B. Ramos-Alvarado, S. Kumar and G. P. Peterson, Phys. Rev. E, 2016, 93, 023101.
156. L. Maffioli, J. P. Ewen, E. R. Smith, S. Varghese, P. J. Daivis, D. Dini and B. D. Todd, Comput. Phys. Commun., 2024, 300, 109205.

---