Regularization of the slip length divergence in water nanoflows by inhomogeneities at the Angstrom scale

We performed non-equilibrium Molecular Dynamics simulations of water flow in nano-channels with the aim of discriminating {\it static} from {\it dynamic} contributions of the solid surface to the slip length of the molecular flow. We show that the regularization of the slip length divergence at high shear rates, formerly attributed to the wall dynamics, is controlled instead by its static properties. Surprisingly, we find that atomic displacements at the Angstrom scale are sufficient to remove the divergence of the slip length and realize the no-slip condition. Since surface thermal fluctuations at room temperature are enough to generate these displacements, we argue that the no-slip condition for water can be achieved also for ideal surfaces, which do not present any surface roughness.


The debate about slip length divergence
In the pioneering paper 1 , it was shown that the Navier slip boundary condition u s = ℓ sγ can be regarded as the low-shear limit of a universal, nonlinear relation between slip velocity u s , the slip length ℓ s , and the local shear rateγ, that presents a divergence of the slip length upon approaching a critical value of the shear rate. The onset of divergence is observed atγ ≃ 10 10 s −1 , both for argon 1 and for water 2 , a value which is accessible so far only to numerical simulations, although the interest in the regime of ultrahigh shear rates is rapidly growing: commercial, industrial-grade fixed-geometry fluid processors like the microfluidizer (Microfluidics Corp, MA, USA) operate at shear rates exceeding 10 7 s −1 to achieve uniform particle size reduction in emulsions and dispersions, to create nanoencapsulations or to produce cell rupture. Rheometers such as the piezoaxial vibrator 3 and the torsion resonator 4 are also employed to measure the viscosity of fluids at rates up to 10 7 s −1 , e.g., for the characterization of ink jet fluid rheology 5,6 . Laser-induced shock-waves in confined water 7 , with speed up to Mach 6 in a 5 µm thick water slab yield even larger numbers, which are of the order of 10 9 s −1 . A more detailed knowledge of the properties of fluids at ultrahigh shear rates has therefore become a compelling task. The issue of slip length divergence, however, is of special interest, as it provides testbed for our understanding of the fluid/solid interactions, and of the minumum requirements for the stability of the fluid slip.
In the approximately fifteen years since the publication of work 1 , violations [8][9][10] , as well as confirmations [11][12][13][14][15] of the ex-  istence of a slip length divergence have been reported, preventing the emergence of a clear consensus on the dependence of fluid slippage properties not only on the imposed shear, but also on the static and dynamic surface inhomogeneities. This lack of consensus calls for further investigations, since controlling and predicting the behavior of water in nanochannels, where surface properties have a profound effect on flow dynamics, stands as a major scientific challenge with many practical applications in modern science and technology. Indeed, an accurate understanding of nanoscale friction phenomena at fluid-solid interfaces is paramount to the design of micro and nanofluidic devices aimed at optimizing mass transport against an overwhelming dissipation barrier [16][17][18][19][20] .
The effects of static properties of the fluid-solid interface, such as the contact angle and surface roughness, on slippage phenomena have been investigated in depth in the recent past 1,10,14,15,[21][22][23][24][25][26][27] , not to mention the growing interest in microtextured and superhydrophobic materials [28][29][30][31][32][33][34][35][36][37] . Recent simulations, for example, suggested that the correlation between hydrophobicity and high slippage could be less strong than expected so far, since some model hydrophilic surfaces can show slip behavior typical of hydrophobic ones 25 . Even though mesoscopic methods have been shown to reproduce quantitatively the results of the molecular dynamics of simple liquids 38 , the description at atomistic level of both fluid and surface is crucial in these kind of investigations, as even very small changes in the chemical or physical properties of the surface can result in remarkably different static and dynamic properties of the fluid, like, for example, in the case of hydroxylated surfaces 39,40 . The role of dynamic properties, such as the wall flexibility and thermal conductivity, remains however much less explored 9,41,42 . Most importantly, the geometric effects due to inhomogeneities of flexible walls, have never been disentangled systematically from the dynamic ones. A recent analysis pointed out that all investigations to date which report absence of a divergent slip length under large shear were performed with flexible walls capable of exchanging momentum and energy with fluid 43 . From this observation, and supported by additional numerical calculations and analytical models, the authors of Ref. 43 concluded that heat and momentum transfer from the fluid to the wall are responsible for the observed saturation of the slip length at high shear rates, as opposed to the divergent behavior which takes place using rigid walls. This picture has been questioned in 44 , where the authors attributed the divergence to the use of a thermostat applied to the fluid molecules. Instead of a saturation of the slip length at high shear rate, these authors report a slippage drop toward the non-slip condition due to fluid heating. With the amount of slippage at high shear rates being reported by different authors to either vanish, reach a constant finite value, or to diverge, the present picture appears rather controversial. More importantly, no clear consensus has yet emerged on the role of different physical quantities governing the transition from a finite to a divergent slip length.
We present here the results of a series of simulations performed with the aim of discriminating those contributions to the slip length which are of purely dynamical origin, from those whose origin is rooted in the surface configurational properties. To this end, we have performed simulations of water flow in channels with flexible walls varying the masses m of the wall atoms and, independently, the harmonic constant k of the potential that bounds wall atoms to their lattice sites (see Fig. 1). According to our findings, this approach proves key to shed new light on the slip length behavior at high shear rates. More precisely, we report numerical evidence that it is the presence of an even extremely tiny amount of disorder in the atomic positions at the surface, rather than wall flexibility, which proves responsible for taming slip length divergence and making it shear-independent.

Non-equilibrium Molecular Dynamics simulations of shear flow
Using a in-house modified version of the Gromacs suite 45 in order to impose Couette flow, we simulated a slab of water molecules confined between two graphite-like walls, each consisting of three atomic layers. The water molecules are modelled using the SPC/E potential 46 , and the electrostatic interaction is calculated using the smooth particle mesh Ewald method 47 , together with Yeh-Berkowitz correction 48 to remove the contribution from periodic images along the wall normal. The distance between the two atomic planes in direct contact with water defines the channel width L = 6.306 nm, while the entire, rectangular simulation cell is 13.035 by 16.188 by 8.0 nm along x, y and z, respectively, with the wall surface normals along the z-direction. The integration time step for the equations of motion was ∆t = 1 fs, and the Couette flow was realized by displacing the positions of the atomic walls every timestep by a fixed amount ±u w ∆t, in addition to the usual movement arising from the integration of the equation of motion. In this way, a constant walls speed ±u w for the upper and lower wall, respectively, is obtained. At the stationary state, the velocity profile (see Fig. 3) is linear over almost the complete extension of the channel, and deviates from it only at distances from the wall which are smaller than the size of a water molecule. By fitting the slope of the fluid velocity, i.e. the effective shear rate,γ eff = ∂ v x /∂ z it is possible to obtain the slip length from its geometrical definition of the distance from the wall at which the (extrapolated) fluid velocity |∂ v x /∂ z|(L/2 + ℓ s ) is equal to the imposed wall speed u w , It has to be noted that the effective shear rate γ eff depends implicitly on the channel width L through the fluid/surface interaction properties, so that Eq. (1), rather than an expression of the slip length in terms of known quantities, is merely a conversion between different ways (slip length and effective shear rate) of measuring the slippage phenomenon: if the slip length is found, e.g., to be independent on the channel width, then the effective shear rate must depend on it so to satisfy Eq. (1), and vice versa. Note also that in this work the symbolγ will be used to denote the imposed shear rate,γ = 2u w /L, not to be confused with the effective shear rateγ eff that develops in the fluid.
The viscosity and friction coefficient can be expressed in terms of the force acting on the surface atoms f as η = f /γ eff and λ = f /u w 20 , respectively, but their ratio is independent from f and equal to u w /γ, which allows to interpret the slip length as a measure of the balance between frictional and viscous forces in analogy with Eq. (1), where the slip length can be expressed through the ratio between imposed and effective shear rate. Likewise, the dependence on the channel width L is not as straightforward as it seems, since the friction coefficient λ depends on the number of water molecules interacting with the wall and, in turn, on the channel width. It could be argued that the appearance of the quantity L in Eqs. (1) and (2) introduces an element of indetermination in the definition of the slip length if the walls are rough, as in this case the channel width (that is, the position of the hydrodynamic boundary) is not well defined. Partial solutions to this problems have been reported in literature. For example, in a work that introduces tunable slip length in DPD simulations 49 , and that has been successfully applied in the simulation of, e.g., nanoconfined electrophoresis 50 and electroosmotic flows 51 , it was pointed out that the slip length and the position of the unknown hydrodynamic boundary can be determined by performing both a Couette and a Poiseuille experiment. In the present case, however, the situation is complicated by the fact that in a Poiseuille flow, in contrast to the Couette flow, the whole spectrum of shear rates is probed, and this would have added a possibly even more detrimental indetermination, since in the present work we focus on the dependence of the slip length on the shear rate itself. By and large, using a fixed position for the hydrodynamic boundaries is probably the best compromise, although in this way the slip length has to be considered as an effective scale that incorporates a possible shift of the real hydrodynamic boundary: higher order corrections could be probably introduced by replacing the geometrical width of the channel with the effective water slab thickness obtained from the density profile. However, these corrections lie outside the scope of this paper.
We investigated two different setups, with respect to the wall-water interatomic interaction potential. In the first setup, hereafter identified as "standard", the interaction between wall and the oxygen atoms in the water molecules is prescribed by a Lennard-Jones potential U(r) = C 12 r −12 − C 6 r −6 with C 6 = 2.47512 × 10 −3 (kJ/nm 6 ) and C 12 = 2.836328 × 10 −6 (kJ/nm 12 ) 52 up to an interatomic distance r of 0.9 nm. Above that distance, the potential is smoothly switched to zero at r = 1.2 nm, using a fourth order polynomial. The second setup is a random quenched functionalization of the first, realized by making 40% of wall atoms purely repulsive (C 6 = 0), and increasing the interaction strength of the remaining ones by a factor α = 1.977. This generates an intrinsic, mainly 1-6 | 3 horizontal inhomogeneity of the surface, e.g. see panel (b) in Fig. 1. The value of α has been chosen such that water equally wets the surfaces of the two setups, resulting in comparable macroscopic contact angles 53 . Given the known dependence of slip length on contact angle 24 , equal wetting is prerequisite to comparing the slip length along two microscopically different surfaces and understanding to what extent the slip length is affected by different types of surface inhomogeneities, rather than by the contact angle.
For each of the surface type just described, we simulated a flexible and a rigid variant. The flexible variant is realized by tethering the atoms of each wall's outermost layer to their respective lattice sites via a harmonic potential U k (r) = 1 2 kr 2 . The spring constant k defines the characteristic excursion ξ = k B T /k of the tethered atoms, so that a small value of k leads to a high effective inhomogeneity. The density profiles of water molecules and surface atoms across the channel for two selected values of ξ are shown in Fig.3. Both profiles are more sharply peaked in the small excursion case, but the effect is significant only next to the surface, and diminishes greatly by the second density peak. In both cases the flow velocity is well represented by the Couette flow solution, even in regions near the wall where the density of water becomes negligible.
The rigid variant, on the contrary, is realized by freezing all surface atoms either at their lattice site position (corresponding effectively to ξ = 0), or in a configuration taken from the equilibrium trajectory of the flexible variant with finite ξ , which we will denote by m = ∞. Both mobile and fixed wall atoms interact with water only, and not among themselves. By varying the mass m of the tethered molecules, it is possible to achieve a separate control over dynamic properties, as the value of m does not influence configurational properties (e.g.: particles distribution, free energies, contact angles), but only dynamical ones, like water-surface friction and thermal conductivity. As a result, we are able to assess the dependence of the slip length on dynamic, static and nonequilibrium properties, controlled by the mass and mean excursion of the wall atoms and by the imposed shear. While this model is surely not a completely realistic representation of a real surface, it allows us to explore in general terms the dependence of slip length as a function of the magnitude of surface inhomogeneities. It is worth noticing that in this work we investigate fluctuations with an extension that is not large enough to model what usually goes under the name of surface roughness. Instead, we are reaching with continuity the limit of very small surface fluctuations, such as those caused by thermal motion. Fig. 1 shows the potential energy isosurfaces for three different cases. While the spatial functionalization (b) already introduces a noticeable perturbation of the energy landscape, the addition of an even tiny spatial displacement of the surface atoms (c) generates the largest perturbation.

The slip length divergence and its regularization
As a first step, with the mass of surface atoms held at a constant m = 1.2 × 10 3 amu (to efficiently integrate the equations of motion of the tethered atoms even with the largest spring constant) we investigated the shear rate dependence of the slip length for different displacements ξ . Providing useful information about the scales which come into play in the determination of slippage properties, we present an overview of this dependence for both the standard and randomized surfaces in Fig. 2. Five different cases were explored for the random quenched surface, between the values ξ = 0 and ξ ≃ 0.11 nm (k = 200 kJ/mol nm −2 ), the former mimicking an infinite spring constant, and the latter corresponding to the one proposed as optimal in 43 . The smallest but still finite value, ξ ≃ 0.0027 nm, was chosen to model a nearly flat wall which yet retains dynamical features. When the wall atoms are perfectly flat (ξ = 0), both surface types show a tendency towards a larger -in fact divergentslip length at increasing shear. It takes only sub-nanoscopic excursions of wall atoms to remove the divergence and make the slip length shear independent. The slip length is shearindependent (constant) also for large values of ξ , however the constant value does increase with the spring constant k (decreasing excursion ξ ) to a limiting value of about 8 nm for the random quenched surface and a slightly larger one for the standard surface. The maximum limiting value for shearindependent slip lengths is reached at ξ c ≃ 0.01 nm. For smaller mean excursions, the slip length is no longer shearindependent and a divergent-like behavior is observed instead. At first glance, a characteristic excursion ξ = 0.01 nm might appear to be surprisingly small, but we note that several other important scales in this system, such as the width of the 1k B T energy basin of a water molecule in the direction normal to the surface (see Fig. 1), and the distance δ ≃ 0.03 nm travelled Lower panel: water viscosity as a function of shear rate. The experimental value η exp = 515.5 kJ/mol/ps/nm 3 at 300K has been interpolated from reference data 55 . by a water molecule during its translational relaxation time of τ r ≃ 0.1 ps 54 . These, however, are only considerations about the orders of magnitude of several process occurring in the system, and a true microscopic explanation of the reason why a length of 0.01 nm is sufficient to trigger the transition from a divergent to a constant slip length is missing.
The transition from shear-independent to divergent slip length is made even more evident by plotting the slip length as a function of ξ (Fig. 4) for two selected shear rates,γ ≃ 0.016 ps −1 at the upper end of the low shear-rates plateau, anḋ γ ≃ 0.08 ps −1 well into the shear region where divergence is observed. At the lower shear rate, saturation is reached for excursions smaller than ξ c ≃ 0.01 nm, whereas the high-shear rate slip length keeps increasing, seemingly unaffected.
From the analysis presented thus far, the role of wall flexibility remains unclear, as this dynamical aspect was not separated out from the other features. For this reason, we selected two extreme values of excursion (ξ = 0.0027 and 0.11 nm) and calculated the slip length for several values of the tethered masses m, including the case of atoms fixed in one equilibrium configuration, m = ∞, as described before.
Results from these calculations are shown in Fig. 5. At low shear rates (γ < 10 −2 ps −1 ) both values of ξ demonstrate marked changes to the slip length with the mass m. The small excursion case ξ = 0.0027 nm halves in slip from the lowest mass m = 1.2 × 10 3 amu to m = ∞, and in the ξ = 0.11 nm case the slip length is reduced to zero at m = ∞. The transition from thermally insulating walls (m = ∞) to conducting walls is therefore associated to a substantial increase in slip length. In other words a flexible wall reduces friction on the fluid (therefore, increasing the slip length) by being capable of adjusting locally to the flux of water molecules. In the extreme case of m = ∞, the slip length divergence is indeed removed, provided that the excursion ξ is larger than ξ c . This provides evidence that the flexibility can not be responsible for the regularization of the slip length at high shear-rates, and acts instead in the opposite direction, by increasing, instead of decreasing, the slip length.
It is important to note that these are entirely surface-related effects. The viscosity of water η is the only bulk property that appears in the definition of the slip length for a Couette flow, Eq. (2), and it does not show any dependence on shear, or spring constant, as shown in Fig. 5, bottom panel. Therefore, we argue that the transition from a divergent to a shearindependent slip length is determined not by the wall flexibility, but rather by the configurational disorder induced by tiny displacements of the wall atoms in the direction normal to the flow.
For Lennard-Jones liquids, Priezjev already noticed 9 that very small, static displacements of the atomic surface atoms (already 7% of the Lennard-Jones diameter) can have a dramatic influence on the slip length using thermal, random or periodic displacements. The systematic analysis performed here puts previous results in perspective, regarding the important case of water slippage. In particular, to get an appreciation for the physical realism of the characteristic excursion lengths discussed here, we compare them to displacements obtained from the effective spring constants of real graphite at 300 K measured by inelastic X-ray scattering 56 , which correspond to about 0.006 nm for nearest neighbors, and are much larger for second neighbors. Thermal fluctuations alone are therefore expected to be sufficient to attain a shear-independent slip length for water.

Conclusions
We have performed non-equilibrium Molecular Dynamics simulations of water flow in nano-channels with separate control of the flexibility and static inhomogeneity of the wall. The simulations show that the disappearance of the divergence of the slip length at high shear rates, formerly ascribed to wall flexibility, is due instead to the excursion of wall atoms from the ideal horizontal plane. Moreover, these results show that atomic displacements as small as those due to thermal motion at room temperature are sufficient to regulate the divergence of the slip length. We conclude that, even in the absence of surface imperfections such as point defects or dislocations, highshear water flows in nanoscale channels should not exhibit any divergent slip length; that thermal motion of the wall atoms is sufficient to tame such divergence.
We acknowledge M. Chinappi & D. Lohse for useful