Slip due to kink propagation at the liquid–solid interface

Abstract

In Couette flow, the liquid atoms adjacent to a solid substrate may have a finite average tangential velocity relative to the substrate. This so-called slip has been observed frequently. However, the particular molecular-level mechanisms that give rise to liquid slip are poorly understood. It is often assumed that liquid slip occurs by surface diffusion whereby atoms independently move from one substrate equilibrium site to another. We show that under certain conditions, liquid slip is due not to singular independent molecular motion, but to localized nonlinear waves that propagate at speeds that are orders of magnitude greater than the slip velocity at the liquid–solid interface. Using non-equilibrium molecular dynamics simulations, we find the properties of these waves and the conditions under which they are to be expected as the main progenitors of slip. We also provide a theoretical guide to the properties of these nonlinear waves by using an augmented Frenkel–Kontorova model. The new understanding provided by our results may lead to enhanced capabilities of the liquid–solid interface, for heat transfer, mixing, and surface-mediated catalysis.

1. Introduction

Slip velocity can be defined as the average tangential velocity of the liquid atoms adjacent to a solid surface. This average motion has been observed frequently (Barrat & Bocquet 1999; Lauga, Brenner & Stone 2007). However, observations of and knowledge about average motion conceal the episodic and short-lived events that constitute slip. We report on transitory groups of atoms with rapid mobility over the solid surface. These localized groups can be the predominant source of slip, rendering the contribution from surface diffusion negligible. In this paper, we investigate and characterize this unheralded mechanism of slip.

We show that at any instant, most liquid atoms adjacent to the solid do not contribute to the slip velocity. Rather, it is only a few groups of liquid atoms that coordinate their motions into fast-moving localized nonlinear waves that convect mass. Under some conditions, the slip velocity is due nearly entirely to the sum of the displacements accrued by these localized waves. In these cases, slip is due not to surface diffusion, but to wave propagation. In their review of slip, Lauga et al. (2007) emphasize the complexity of phenomena that take place at the liquid–solid interface. With the results presented here, we add another element to their array of interesting phenomena that contribute to slip.

There have been many elegant and artful numerical, analytical and experimental studies of liquid slip, summarized in § 2, including molecular dynamic simulations similar to our simulations, which are described in § 3. We find a previously unobserved phenomenon of rapidly propagating waves at the liquid–solid interface. This finding is made possible by the insight that the solitonic nature of localized disturbances produces a unique signature. That signature is the doubly occupied lattice cell, made visible by measuring the locations of the liquid atoms relative to the substrate lattice, as we describe in § 4. We can then scan our data for the occurrence of a doubly occupied cell in order to extract those events most meaningful for slip. In particular, the doubly occupied cell is central to a localized nonlinear wave that transports mass and leads to slip, as described in § 5. Observations from our molecular dynamics (MD) simulations agree well with the predictions of the Frenkel–Kontorova and sine-Gordon equations, lending analytical support characterizing our observations as nonlinear waves akin to solitons, as further discussed in § 5. We summarize the main points of our study and note their applicability in the final § 6.

2. Prior research

The no-slip boundary condition remains a dependable boundary condition for many macroscopic Newtonian flows. However, when liquids are confined to channels of a width of only a few molecular diameters, the conventional hydrodynamic theories based on the Navier–Stokes equations may fail, and the no-slip boundary condition may no longer be valid (Travis, Todd & Evans 1997; Travis & Gubbins 2000; Zhu & Granick 2002). Experimental, analytical and numerical studies showed the now well-known slip: that liquids adjacent to a solid move with a finite velocity relative to the solid surface (Lauga et al. 2007; Shu, Teo & Chan 2018; Wang et al. 2021). Recently, Hadjiconstantinou (2024) presented a discussion of the state of the art of the theoretical understanding of slip. Here, we further discuss prior works that are most relevant to our study.

Together with the equations that now bear his name, (Navier 1822) derived what became known as the Navier slip condition:

(2.1)\begin{equation} E U + \epsilon\,\frac{{\rm d} U}{{\rm d} z} = 0, \end{equation}

shown here for the tangential component of the liquid velocity $U$ above a flat wall whose normal is in the $z$-direction. Navier delegated constants $E$ and $\epsilon$ as sliding resistances: $\epsilon$ as the sliding resistance between fluid molecules in adjacent layers, the present-day dynamic viscosity, and $E$ as the sliding resistance between the molecules in the vicinal liquid layer and the wall itself, the present-day liquid–solid friction coefficient. The ratio $\epsilon /E$ is the slip length. Navier's presentation was in terms of molecules moving in layers, a point of view that we also use in defining the first liquid layer (FLL) in § 3. In fact, near a wall, the liquid density is layered (Rowley, Nicholson & Parsonage 1976; Abraham 1978). Careful formulations and analyses of the slip boundary condition can be found in Miksis & Davis (1994), Brenner & Ganesan (2000) and Bolaños & Vernescu (2017), including over textured and patterned surfaces (Priezjev & Troian 2006; Kamrin, Bazant & Stone 2010; Luchini 2013; Zampogna, Magnaudet & Bottaro 2019).

Koplik, Banavar & Willemsen (1989) provided an early yet thorough MD simulation of liquid slip, including examination of Couette and Poiseuille flows, contact line motion, and the trajectories of the individual liquid atoms. Our MD simulation methodology is similar to theirs in that we make use of molecular trajectories to identify the atomic-level mechanisms of liquid slip, as described in § 4.

Models of surface diffusion at the liquid–solid interface have presented the dependence of slip on shear rate, temperature and molecular interaction parameters (Hsu & Patankar 2010; Shu et al. 2018; Wang & Hadjiconstantinou 2019; Wang et al. 2021; Shan et al. 2022). Hadjiconstantinou (2021) formulated and solved the Fokker–Planck dynamics of a point particle moving in one dimension, driven by shear, slowed by friction with a periodic substrate, and interacting with a thermal background. As shear force is increased, these results reveal slip length at first insensitive to shear, then rapidly increasing to a higher value that once again is invariant to shear. As pointed out by Hadjiconstantinou, these general features of the slip length are consistent with those that were found by other means (Martini et al. 2008b).

Motions more complex than simple diffusion have been anticipated to contribute to the surface transport at the liquid–solid interface. In particular, several types of mechanisms have been described, whereby atoms’ correlated motion with neighbouring atoms has been described as gliding, reptation and dislocation motions (Oura et al. 2013). Studies of monolayers range from numerical studies of atoms on atomically structured substrates (Çam, Lichter & Goedde 2021) to micron-sized beads over a laser-generated interference pattern (Reguzzoni et al. 2010; Brazda et al. 2018). These studies reveal that even at low levels of forcing, patches of molecules will slip due to the incommensurate alignment between the particles and the surface. In contrast, observations of coordinated motion at the liquid–solid interface under macroscopic flows such as Couette or Poiseuille flow are lacking. It is precisely to seek such coordinated motion that we undertook this investigation.

3. The MD simulations of liquid slip

Our non-equilibrium MD simulations are performed using the nanoscale MD (NAMD) programme (Phillips et al. 2005). We simulate a liquid confined between two thermal solid walls, shown in figure 1. Each wall consists of six $20\times 20$ layers of atoms in the $x$- and $y$-directions. The solid atoms are maintained on a cubic lattice with lattice spacing $\lambda$ by nearest-neighbour linear springs of stiffness 20 N m$^{-1}$. Periodic boundary conditions are implemented along the $x$- and $y$-directions of the square domain, which has dimensions $L_{x} = L_{y} = 20\lambda$. Periodic boundary conditions are also used in the $z$-direction, with the periodic box size chosen to be large enough that the periodic images of the upper and lower walls are always separated by a distance greater than the NAMD non-bonded interaction cutoff distance 1.4 nm. The channel width $h$, equal to the $z$ distance between the mean positions of the innermost solid layers, is 18 liquid molecular diameters.

Figure 1. Snapshot of the Couette flow set-up. Liquid and solid atoms are shown as cyan and yellow spheres, respectively. Solid yellow lines indicate the linear bonds between the nearest-neighbour solid atoms. The direction of motion of each solid wall is shown by a large yellow arrow. The channel width $h$ is the $z$ distance between the mean positions of the innermost solid layers.

The non-bonded interaction potentials are given by

(3.1)\begin{equation} \phi_{ij}(r) = 4\epsilon_{ij} [(\sigma_{ij}/r)^{12} -(\sigma_{ij}/r)^{6}], \end{equation}

where $\phi _{ij}$ indicates the Lennard-Jones interaction strength between atoms of types $i$ and $j$. Indices $i$ and $j$ are chosen to indicate the atomic species, liquid (L) or solid (S). The Lennard-Jones potential provides a simple model of atomic interaction, with only an energy parameter $\epsilon _{ij}$ and a length parameter $\sigma _{ij}$. The energy and length parameters associated with the interactions between the liquid and solid atoms are determined using the Lorentz–Berthelot combining rules (Allen 2004):

(3.2a,b)\begin{equation} \sigma_{LS} = \frac{\sigma_{LL} + \sigma_{SS}}{2} \quad \mathrm{and} \quad \epsilon_{LS} = \sqrt{\epsilon_{LL}\epsilon_{SS}}. \end{equation}

Table 1 summarizes the energy and length parameters, and the masses of the liquid and solid atoms. In all simulations, the length parameter for the solid–solid interaction is chosen such that the non-bonded equilibrium spacing of the solid atoms, $r_{0}^{SS}\equiv 2^{1/6}\sigma _{SS}$, is equal to the solid lattice spacing $\lambda$.

Table 1. Lennard-Jones energy $\epsilon _{ij}$ and length $\sigma _{ij}$ parameters associated with the liquid–liquid (LL), solid–solid (SS) and liquid–solid (LS) interactions as well as the masses of liquid and solid atoms. Here, $k_{B}$ is the Boltzmann constant.

In some MD simulations, the solid atoms bounding the liquid flow are fixed in their lattice positions and not given the freedom to fluctuate, constituting a so-called rigid or athermal wall. Martini et al. (2008a) have shown that at high shear rates, slip is more accurately reproduced by allowing the wall atoms to possess their thermal energy and fluctuate, as we do in this study. Furthermore, liquid atomic trajectories are most faithfully simulated when the viscous heat generated by the sheared liquid is allowed to flow to the solid boundaries and, only there, be thermostatted (Bernardi, Todd & Searles 2010; Yong & Zhang 2013). Therefore, in this study, we thermostat neither the liquid itself nor the adjacent solid layers, but thermostat only the outer solid layers most distant from the liquid.

The simulation protocol consists of three stages. In the first stage, the walls are not yet set into motion while the entire system is coupled to the Langevin thermostat. This first stage allows the liquid atoms to diffuse into random positions. In the second stage, the thermostat is restricted to the two outer solid layers in each wall. During this stage, the top and bottom solid walls accelerate to their target velocities $U_{WALL}$ in the $x$-direction, as shown by the horizontal yellow arrows in figure 1. Target velocities of the top and bottom solid walls are equal in magnitude but opposite in direction. In the third stage, both solid walls move with constant velocity in the $x$-direction. We report only the results obtained during the third stage, in which the average properties of the liquid and solid remain constant.

As discussed in detail below, we carry out two parametric studies. In one, the wall velocity $U_{WALL}$ is systematically varied from 0.2 to 220 m s$^{-1}$. For these studies, we present the slip velocity as a function of shear rate. We calculate the shear rate by finding the slope of a straight-line fit to the liquid velocity profile in $z$-direction near the centre of the channel.

In the second study, we present results where the Lennard-Jones size of the liquid atoms $r_{0}^{LL}\equiv 2^{1/6}\sigma _{LL}$ is kept constant while the lattice spacing $\lambda$ and the solid–solid length parameter $\sigma _{SS}$ are varied, yielding the ratio $r_{0}^{LL}/\lambda$ of the relative size of the liquid atoms to the solid atoms as the control parameter.

For data collection at sampling rate 1 ps, the duration of the third stage is 4–26 ns, depending on the shear rate. For each value of $r_{0}^{LL}/\lambda$, we run four or five simulations with distinct initial positions and velocities for the liquid atoms, while all other simulation parameters remained the same. Sampling at 1 ps would be sufficient if our only objective is to measure the average properties of the system such as temperature, density, slip velocity and slip length. However, it typically takes less than a picosecond for a liquid atom to move from one substrate equilibrium site to the next. Consequently, for our goal to resolve the detailed trajectories of the coordinated motion of individual atoms, femtosecond resolution in atomic positions is needed. In simulations where we gather data on these dynamics, we use sampling rate 10 fs and run the simulations up to 2 ns.

A well-established observation of atomic-level structure is the variation in liquid density near a solid (Rowley et al. 1976; Abraham 1978; Barrat & Bocquet 1999; Morciano et al. 2017). We too find this layering, and make use of it in our analysis. In particular, we define the FLL as the volume between the mean positions of the innermost solid layer and the first minima of the liquid density profile, as shown in figure 2. Liquid atoms whose positions lie within this volume constitute the FLL. Far enough from the solid walls, the variations disappear and the liquid density becomes constant. The total number of liquid atoms in the channel is chosen such that the liquid density, denoted by $\rho$, has the constant value $\rho =0.87\sigma _{LL}^{-3}$ at the centre of the channel.

Figure 2. The liquid density profiles $\rho \sigma _{LL}^{3}$ as functions of the $z$-position of the liquid atoms $z/\sigma _{LL}$. Here, $z$ is measured relative to the innermost solid layer. Blue, green and red solid lines indicate the liquid density profiles of $r_{0}^{LL}/\lambda = 1.54$, 1.00 and 0.84, respectively. The first minimum defining the upper boundary of the FLL is nearly the same for all densities, and is shown by a dotted blue line for $r_{0}^{LL}/\lambda =1.54$. Far enough from the solid walls, the liquid densities asymptote to a common value $\rho \sigma _{LL}^{3}=0.87$.

Averaging the velocities relative to the moving substrate of all FLL atoms in the $x$-direction over a long enough interval of time yields the slip velocity $U_{FLL}$, shown in figure 3. At lower shear rates, longer simulations are required to resolve the ever-decreasing slip velocity. Simulations at zero shear rate yield zero slip velocity approximately ${\pm }$1 m s$^{-1}$. Thus from zero shear rate to all but the largest shear rates shown in the inset of figure 3, our data are consistent with the proportionality of shear rate and slip velocity. The slopes of the straight-line fit to the data shown in figure 3 yield the Navier slip lengths $(6.47, 0.54, 0.84) \sigma _{LL}$ for the three values of $r_0^{LL}/\lambda$, $(1.54, 1.0, 0.84)$, respectively.

Figure 3. The slip velocity $U_{FLL}$ of the FLL, as a function of the shear rate $\dot {\gamma }$. Blue circles, green triangles and red squares indicate the slip velocities as the substrate wavelength $\lambda$ increases relative to the size $r_0^{LL}$ of the liquid atoms, i.e. $r_{0}^{LL}/\lambda = 1.54$, 1.00 and 0.84, respectively. The dashed lines are linear fits to $U_{FLL}$ for shear rates $\dot {\gamma }<10^{10}\,\mathrm {s}^{-1}$. The slopes of these lines yield Navier slip lengths $(6.47, 0.54, 0.84) \sigma _{LL}$ for blue, green and red, respectively. The inset shows $U_{FLL}$ over a larger range of shear rates.

As also shown in the figure, the slip length depends on the lattice spacing of the substrate, with the slip length increasing as the lattice spacing of the substrate decreases. Other researchers have shown the dependence of slip velocity (and slip length) on the geometrical and chemical properties of the substrate (Thompson & Troian 1997; Barrat & Bocquet 1999). In general, then, we find that slip velocity is not a universal value independent of the spacing of the atoms along the surface, but is dependent on the atomic type and the wavelength of the solid substrate.

Figure 4 shows the dependence of slip velocity on $r_0^{LL}/\lambda$ for wall speed $U_{WALL}=180$ m s$^{-1}$. Ratios $r_0^{LL}/\lambda <1$ indicate that the substrate lattice spacing is greater than the equilibrium spacing of the liquid atoms. Conversely, ratios $r_0^{LL}/\lambda >1$ indicate that the equilibrium spacing of the solid atoms is smaller than the equilibrium spacing of the liquid atoms. Each data point represents the average of five simulations of duration 6 ns; the error bars represent the standard deviations of each set of five simulations. A minimum in the slip velocity occurs at $r_{0}^{LL}/\lambda =1$, where the substrate wavelength is equal to the equilibrium spacing of the liquid atoms.

Figure 4. A close-up of the slip velocity $U_{FLL}$ of the FLL for $r_{0}^{LL}/\lambda$ near unity. The slip velocity has a minimum at $r_{0}^{LL}/\lambda =1$. Error bars are the standard deviations of five simulations for each $r_{0}^{LL}/\lambda$. Atomic positions and velocities are sampled every 1 ps. Wall speed is $U_{WALL}=180$ m s$^{-1}$. The inset shows $U_{FLL}$ across a wider range, $0.78 \leq r_{0}^{LL}/\lambda \leq 2.22$.

4. Mechanisms of liquid slip

The simple cubic arrangement of the substrate atoms creates a square-symmetric potential energy landscape through which the liquid atoms move, as illustrated in figure 5. On a plane at a small height above the substrate, the largest values of the potential energy are found directly above the substrate atoms. Conversely, the minimum values of the potential energy are located equidistant from four neighbouring solid atoms. We denote these locations, shown as white crosses in the figure, as substrate equilibrium sites. Due to this pattern in the potential energy landscape, the atoms in the FLL are most likely to be found near the substrate equilibrium sites, and they are least likely to be found directly over the locations of the solid atoms. As a result, it is unusual to find more than one liquid atom crowded into the neighbourhood of an equilibrium site. Indeed, a doubly occupied neighbourhood of an equilibrium site is the critical event for localized slip, as we will show.

Figure 5. Potential energy landscape as a function of the $x$- and $y$-positions of a liquid atom at a fixed height above the substrate. The strength of the potential energy is colour-coded as dark blue for low and bright yellow for high potential energy locations. White cross markers indicate the substrate equilibrium sites. Horizontal and vertical white dashed lines indicate low-energy corridors connecting neighbouring equilibrium sites. Filled black circles indicate the lattice positions of the solid substrate atoms. Solid vertical and horizontal black lines indicate the boundaries of substrate cells, as defined in the text. A patch of only 9 cells of 400 is shown here.

Using the pattern shown in figure 5, we divide the FLL into 400 cuboidal substrate cells, each centred over a substrate equilibrium site. The base of each cell is bounded by four substrate atoms (the black circles in figure 5) and the (black) lines connecting them. The height of the cells is equal to the height of the FLL. In figure 5, we show a top view of a patch of 9 neighbouring cells (out of 400 total for our $20\times 20$ periodic lattice) superimposed on the substrate potential energy landscape.

Using the data from simulations sampled every 10 fs, we label and identify all the liquid atoms in the FLL, and track the position of each uniquely labelled atom in time, paying particular attention to the cell that each liquid atom occupies. This allows us to follow each atom in the FLL as it moves between substrate cells. Figure 6 shows a schematic view of a typical event as an atom in the FLL hops from one cell to the next. In figure 6(a), each cell is occupied by a single liquid atom, which is a common configuration of atomic positions in the FLL. However, a liquid atom can spontaneously move from its current cell into an already-occupied cell, creating a doubly occupied cell. Figure 6 shows an example in which a liquid atom moves from the cell labelled A to the cell labelled B in the time between Figures 6(a,b), creating doubly occupied cell B, and leaving behind vacant cell A.

Figure 6. Schematic representation of liquid atoms (blue circles) on a small square grid of cells at two consecutive instants of time, shown in (a,b), respectively. Due to the $x$-motion of a liquid atom, shown by the horizontal blue arrow from cell A to cell B, cell A becomes vacant, and cell B becomes doubly occupied. In (a), each cell is occupied by a single liquid atom. In (b), cell A is not occupied by any liquid atom and is a vacant cell, and cell B, occupied by two liquid atoms (each coloured red), is designated as a doubly occupied cell. (The diameter of the liquid atoms is chosen for clarity and does not represent their characteristic size.)

After an atom moves into an already-occupied cell, the original occupant of the now doubly occupied cell can become dislodged and move forward, so the double occupancy moves onwards into the adjacent cell. This process can repeat itself, as shown schematically in figure 7. As the doubly occupied cell propagates, atoms advance by one cell forwards. Therefore, the mean motion of the entire layer can be produced by the few atoms moving quickly and sequentially from one doubly occupied cell to the next, in a coordinated fashion.

Figure 7. Schematic representation of a doubly occupied cell moving in the $x$-direction. A horizontal strip of seven neighbouring cells in the $x$-direction is shown, sampled at four consecutive instants, from (a) to (d). As the doubly occupied cell moves, atoms advance by one cell in the $x$-direction, and the pair of atoms in the doubly occupied cell (coloured red) evolves sequentially. Horizontal red arrows indicate the atoms that advance into the next cell. Also shown by the horizontal blue arrow, in the time between (c) and (d), is the motion of a liquid atom from a singly occupied cell into a vacant cell, resulting in the vacant cell moving one cell to the left.

The motion of a liquid atom into a vacant (unoccupied) cell is quite different from the motion of an atom into an already-occupied cell. When a liquid atom moves into a vacant cell, nearby liquid atoms are nearly unaffected. Therefore, such motion represents an uncoordinated, independent atomic motion over the substrate, which is the typical description of surface diffusion as a mechanism of slip.

We identify and distinguish every instance of an atomic hop either into an already-occupied cell or into a vacant cell. We track the motion in both the $x$- and $y$-directions. However, as there is no net motion in the transverse $y$-direction, we restrict our analysis here to the $x$-direction only.

Using the data from simulations sampled at 10 fs, we separately calculate the net liquid displacement in the $x$-direction due to (i) atomic motion involving doubly occupied cells, and (ii) atomic motion into vacant cells. When divided by the duration of the simulation and the time-averaged number of atoms in the FLL, these displacements provide the fractional contributions to slip due to atomic motion into doubly occupied or vacant cells.

In figure 8, we show the slip velocity of the FLL over the range $0.78 \leq r_0^{LL}/\lambda \leq 2.22$ for wall speed $U_{WALL}=180\,\mathrm {m}\,\mathrm {s}^{-1}$. For each value of $r_0^{LL}/\lambda$, we run a single simulation with sampling rate 10 fs and simulation duration 2 ns. The blue circles show the slip velocity for the layer as calculated directly from the atomic velocities. (This is the same range of data as shown in figure 4.) The red squares show the slip velocity due only to motion involving doubly occupied cells. The cyan diamonds show only the slip due to atomic motion into vacant cells. Taken together, the sum of these two contributions accounts for the totality of the slip. Liquid slip due to doubly occupied cells becomes increasingly predominant as $r_{0}^{LL}/\lambda$ decreases. The cross-over value of $r_{0}^{LL}/\lambda$ at which the contributions to slip are equally divided between doubly occupied cells and atoms hopping into vacant cells occurs near $r_{0}^{LL}/\lambda \simeq 1$. Slip is overwhelmingly due to doubly-occupied cells for $r_{0}^{LL}/\lambda < 1$, and continues to contribute significantly to slip even for values of $r_{0}^{LL}/\lambda$ above unity.

Figure 8. The slip velocity $U_{FLL}$ of the FLL, as a function of $r_{0}^{LL}/\lambda$, at the wall speed $U_{WALL}=180$ m s$^{-1}$. Blue circles indicate the slip velocity as determined by the direct measurement of the atomic velocities. Red squares indicate the slip velocity due to the mass propagation involving doubly occupied cells alone. Cyan diamonds indicate the slip velocity due to the mass propagation by atoms hopping into vacant cells. Atomic positions and velocities are sampled every $10$ fs.

To provide a sense of the role played by vacant and doubly occupied cells in the generation of slip, figure 9 shows an annotated view of the entire FLL obtained from the simulation with $r_{0}^{LL}/\lambda =1$ at wall speed $U_{WALL}=20$ m s$^{-1}$ and sampling rate 10 fs. Approximately 5 % of the liquid atoms are in a doubly occupied cell, while approximately 95 % of the atoms merely fluctuate around their equilibrium sites. Also shown are the trajectories of doubly occupied cells as they were tracked moving over several cells; for these parameter values, the lifetime of the motion involving a doubly occupied cell is typically a few picoseconds.

Figure 9. Portrait of the FLL constructed from the MD data. The blue translucent circles indicate the instantaneous positions of the liquid atoms. Liquid atoms in doubly occupied cells are coloured red. The thick green lines indicate the trajectories over time of doubly occupied cells. The open green squares indicate vacant cells. The black dots indicate the lattice positions of the solid substrate atoms. The mean slip direction is to the right. Here, the wall speed is $U_{WALL}=20$ m s$^{-1}$, $r_{0}^{LL}/\lambda =1$, and the atomic positions are sampled every 10 fs.

To highlight the sequential nature of the atomic motion involving doubly occupied cells, we show in figure 10 the $x$-positions of six neighbouring liquid atoms as a function of time as the atoms hop in the $x$-direction between doubly-occupied sites. The atomic positions shown are obtained from the same simulation that is used to generate figure 9. The actual atomic trajectories are shown in figure 10(a); a schematic representation of the motion is shown in figure 10(b) for clarity.

Figure 10. Atomic trajectories illustrating coordinated atomic motion into or out of doubly occupied cells. (a) Excerpt from our MD data, simplified in (b) to highlight the advance of each liquid atom 1–5 from one equilibrium site to the downstream site. Each atom 1–5 is advanced one substrate lattice spacing in the $x$-direction by the passage of the doubly occupied cell, instantaneously located at the thickened red portion of the atomic trajectories. The thin vertical orange lines mark the first and final times of the occurrence of the doubly occupied cell. The $x$-boundaries between cells are marked by dashed grey lines. Note that the liquid atoms tend to remain near the substrate minima located halfway between adjacent dashed lines. Magenta portions of the atomic trajectories mark when the atom has drifted out of the FLL. Here, the wall speed is $U_{WALL}=20$ m s$^{-1}$, $r_{0}^{LL}/\lambda =1$, and the atomic positions are sampled at every 10 fs.

In this example, atom 1 starts in cell A and hops to cell B near $t=5$ ps, creating a doubly occupied cell, which is indicated by highlighting the trajectories for atoms 1 and 2 in red. The interaction of atoms 1 and 2 results in atom 2 hopping from cell B to cell C (already occupied by atom 3). In this manner, the doubly occupied site propagates from cell B to cell F. The sequence terminates when atom 6 moves into the unoccupied cell G at the top of the diagram. The sequence of doubly occupied cells exists during the time interval between the two vertical orange lines, for approximately 0.5 ps. During this period, the doubly occupied cell translates by five substrate wavelengths, and five atoms propagate by one lattice spacing forwards.

5. Frenkel–Kontorova model of slip

The two atoms in doubly occupied cells are part of a larger number of atoms whose motion is correlated. When properly visualized, this group of atoms is seen to take on a shape that is characteristic of nonlinear waves in the sine-Gordon equation and Frenkel–Kontorova (FK) model (Braun & Kivshar 1998).

We define an offset for each atom near the doubly occupied cell,

(5.1)\begin{equation} \textrm{[offset of atom $k$]}=\frac{x_k(t)-k\lambda}{\lambda}, \end{equation}

where $x_k(t)$ is the $x$-coordinate of atom $k$, and $k\lambda$ is the position of the atom's initial substrate equilibrium site.

Figure 11 shows the offsets of the atoms in the neighbourhood of a propagating doubly occupied cell. Downstream of the doubly occupied cell, atoms gradually move away from their equilibrium sites as they are swept into the doubly occupied cell. As the doubly occupied cell propagates onwards, upstream atoms relax into new equilibrium sites after being transported one cell forward. These profiles, such as shown in figure 11, thereby refer to the distribution of atomic offsets relative to the atoms’ initial equilibrium positions. So offset 0 shows that atom $k$ has not moved away from its initial site, while offset 1 indicates that the atom has been transported forwards by the distance of one substrate cell. Thus the doubly occupied cell marks a larger propagating structure of displaced atoms that we call a kink.

Figure 11. Time-averaged kink profiles derived from atomic positions in the FLL of the MD simulations. Data points of blue triangles, green circles and red squares are the mean values of atomic offsets from MD simulations with $r_{0}^{LL}/\lambda$ for 1.08, 1.00 and 0.87, respectively. The wave profiles are well described by (5.2), shown as solid curves, with steepnesses $A=0.27$ (blue triangles), $A=0.44$ (green circles) and $A=0.85$ (red squares). The location of the doubly occupied cell that marks the centre of the propagating kink is shown by the diamond at the top of the figure.

The profiles seen in figure 11 can be fitted to the sigmoid function

(5.2)\begin{equation} \textrm{[offset of atom $k$]}=\frac{2}{\rm \pi}\tan^{{-}1} \{\exp[{-}A(k-k_0)]\}, \end{equation}

where the parameters $A$ and $k_0$ describe the steepness of the profile and the location of its centre, respectively. This function is the discrete version of the soliton solution to the sine-Gordon equation, which is the continuum limit of the standard FK model (Braun & Kivshar 1998).

For each of the three curves shown in figure 11, the data points are the time-averaged atomic positions over the lifetime of one kink. All frames are centred by aligning the location of the doubly occupied cell. The curves are then generated by fitting the averaged atomic positions to the sigmoidal (5.2) to find $A$ and $k_0$ for each profile.

To understand the shape and behaviour of the kinks observed in the MD simulations, we model the FLL using a modified version of the one-dimensional FK model. The FK model describes $N$ particles along the $x$-axis at positions $x_k(t)$, $k=1, 2, \ldots, N$, that are interacting through a linear spring force with their neighbouring particles, and are also interacting with a sinusoidal force that arises from the periodic potential of the adjacent substrate (Braun & Kivshar 1998). If the spacing between the particles is precisely equal to the wavelength of the potential, then the particles achieve their lowest energy state with each particle at the periodically located potential minima. When an additional particle is present, such that there are one too few minima to host each of the particles, the FK model predicts that the spacing between particles is no longer uniform. While most particles remain periodically placed close to substrate potential minima, a few particles within a narrow spatial interval are crowded together. Thus the FK model can be used to understand the profiles shown in figure 11.

For our application, the particles of the FK model are the liquid atoms, and the potential is due to the substrate of solid atoms. Prior work showed that the FK model provides a description of the dependence of slip on shear rate (Lichter, Roxin & Mandre 2004; Martini et al. 2008b). Furthermore, predictions of the FK model agree with the observations that slip is bounded in the limit of large shear rates over thermalized solid substrates (Martini et al. 2008a). Though these studies gave evidence that the FK model could describe certain aspects of mean slip, there was little comparison with the observations at the molecular scale. This was due to the lack of data on molecular trajectories and the absence of a conceptual framework that permitted the localized propagating structures to be identified and extracted from the noisy MD data. These deficiencies have been overcome in the present work through the recognition of the doubly occupied cell.

To more realistically model the interactions between neighbouring atoms in a liquid, we replace the linear nearest-neighbour coupling of the classical FK equation with a nearest-neighbour force arising from the 12-6 Lennard-Jones potential. In addition, the liquid atoms are forced into motion by a shear force of the form $\eta _{LL}(V-\dot {x}_k)$, where $V$ is the constant velocity of the liquid layer above, and $\eta _{LL}$ is the friction coefficient between the adjacent liquid layers. As the liquid atoms move with speed $\dot {x}_k$ over the substrate, they are also subject to a frictional force proportional to their speed $\eta _{LS}\dot {x}$, where $\eta _{LS}$ is the liquid–solid friction coefficient.

With these modifications, the FK model becomes

(5.3)\begin{align} \ddot{x}_k &={-}\sin{(x_k)} \nonumber\\ &\quad + g\left[\left(\frac{\sigma}{x_{k+1}-x_k}\right)^7- 2\left(\frac{\sigma}{x_{k+1}-x_k}\right)^{13} -\left(\frac{\sigma}{x_k-x_{k-1}}\right)^7+ 2\left(\frac{\sigma}{x_k-x_{k-1}}\right)^{13}\right] \nonumber\\ &\quad +\eta_{LL}(V-\dot{x}_k) - \eta_{LS}\dot{x}_k, \end{align}

where

(5.4a,b)\begin{equation} g=\frac{12\epsilon_{LL}\lambda}{{\rm \pi} \sigma_{LL} \phi_{0}} \quad \mathrm{and} \quad \sigma = \frac{2{\rm \pi}\sigma_{LL}}{\lambda}, \end{equation}

and where the particle positions and time have been non-dimensionalized in the usual way (Braun & Kivshar 1998). The non-dimensional parameter $g$ signifies the strength of the liquid–liquid interaction relative to that of the liquid–substrate interaction. The ratio of the characteristic size of the liquid atoms relative to the lattice spacing, $\sigma$, plays an important role in setting the ground state and the dynamics of the solutions to the FK equation (Aubry 1983; Braun & Kivshar 1998).

The dimensionless parameters defined in (5.4a,b) depend on the dimensional liquid–liquid Lennard-Jones energy $\epsilon _{LL}$, the length parameter of the liquid–liquid interaction $\sigma _{LL}$, the substrate wavelength $\lambda$, and the dimensional amplitude of the sinusoidal substrate potential energy $\phi _{0}$. While the first three can be taken directly from the input parameters for the MD simulations, the substrate potential amplitude used in the FK model must be estimated from the MD interaction potential between all the substrate atoms and an atom in the FLL. We do this by calculating the potential energy above two $(x,y)$-locations in the MD simulations: (i) a substrate equilibrium site, yielding $\phi _{eqb}(z)$; and (ii) halfway between two neighbouring equilibrium sites, yielding $\phi _{jump}(z)$. Atoms in the FLL of the MD simulations are expected to tend towards heights at which the potentials are at a minimum. We therefore take the substrate potential amplitude for the modified FK model to be one-half of the difference between the minimum values of $\phi _{eqb}(z)$ and $\phi _{jump}(z)$, yielding

(5.5)\begin{equation} \phi_{0} = \frac{\min{[\phi_{jump}(z)]} - \min{[\phi_{eqb}(z)]}}{2}. \end{equation}

This results in a substrate potential amplitude for the modified FK model $\phi _{0}/k_{B} \simeq 35.5(1.1-r_{0}^{LL}/\lambda ) + 66$ K.

We use the modified FK model to describe the dynamics of the atoms near a doubly occupied substrate cell. As in our MD simulations, we use periodic boundary conditions in the modified FK model, with periodic domain length $(N-1)\lambda$, where $N$ is the number of particles. This assures that there is an extra particle in the FK chain relative to the number of substrate potential minima. Figure 12(a) shows numerical solutions of the FK equations (5.3) at multiple instants of time. We call these nonlinear wave solutions to the FK model ‘FK kinks’. FK kinks propagate at a constant speed with an invariant profile. The FK kink profiles are closely approximated by the sine-Gordon soliton solution (dotted curves) from the sigmoid function (5.2). The profile can be characterized succinctly by its steepness parameter $A$, which indicates the narrowness of the FK kink profiles, three of which are shown in figure 12(b).

Figure 12. Numerical results from the FK model. (a) The propagation of a single FK kink through 52 contiguous atoms is shown at intervals of 0.7 ps. The two atoms instantaneously inside a doubly occupied cell are highlighted in red, and the locations of their centres of mass are shown by diamonds at the top of the figure. The dotted lines at the initial and final frames are the nearly identical sine-Gordon soliton profiles, fitted with $A=0.27$ for $r_{0}^{LL}/\lambda =1$. (b) Time-averaged FK kink profiles for three different values of $r_{0}^{LL}/\lambda$, showing the steepening profile as $r_{0}^{LL}/\lambda$ decreases: blue triangles for $r_{0}^{LL}/\lambda = 1.08$, $A=0.17$; green circles for $r_{0}^{LL}/\lambda = 1.00$, $A=0.27$; red squares for $r_{0}^{LL}/\lambda = 0.87$, $A=0.56$. Data points are from numerical solutions to the FK model, and correspondingly coloured curves are best fits from the sigmoid function (5.2).

To compare the wave properties of MD kinks to FK kinks, we compare their velocities and narrowness. Figure 13(a) shows the velocity of the kinks, $v_{k}$, in both the MD and FK simulations as a function of $r_{0}^{LL}/\lambda$. We calculate kink velocities in the MD simulations by dividing the kink displacements in the $x$-direction by their lifetimes. It is observed that the kink velocities in MD simulations increase as $r_{0}^{LL}/\lambda$ increases. This trend is also observed in our FK simulations. The agreement between the kink velocities from MD and FK simulations gets better as $r_{0}^{LL}/\lambda$ decreases, where we can sample a larger number of propagating kinks, which is signified by smaller error bars. Figure 13(b) shows the steepness parameter $A$ of the kinks of both the MD and FK simulations as a function of $r_{0}^{LL}/\lambda$. Kinks get narrower as $r_{0}^{LL}/\lambda$ decreases.

Figure 13. Results from the MD simulations compared with those from the FK model. (a) The kink velocity $v_{k}$ as a function of $r_{0}^{LL}/\lambda$. Filled red circles indicate the average MD kink velocities. The thick green line indicates the velocities of the FK kinks obtained from the FK simulations. (b) The steepness parameter $A$ of the kinks as a function of $r_{0}^{LL}/\lambda$. The MD data (filled red circles) were averaged over all kinks that propagated at least four substrate wavelengths. Error bars are the standard deviations of all those kinks for each $r_{0}^{LL}/\lambda$. For the green line of the FK model, we set $V=1$, and use values $\epsilon _{LL}/k_{B}=188\,\mathrm {K}$, $\sigma _{LL}=(2{\rm \pi} )2^{-1/6}$, $\phi _{0}/{k}_{B}\simeq 70\,\mathrm {K}$, identical to those of the MD simulation. The friction factors were chosen to comply with the high shear rate limit of the FK equation, $\eta _{LL}/\eta _{LS}={O}(1)$; see Appendix A.

When compared with the predictions of the modified FK model, both the profile shape, as measured by $A$, and the velocity of the kinks seen in the MD simulations, show similar trends and comparable values as a function of $r_{0}^{LL}/\lambda$. In assessing the agreement between the two sets of data, the extreme simplicity of the FK model should be recalled. In the FK model, particles move in only one dimension and the system is closed. This is to be contrasted with the three-dimensional motion of atoms in the MD simulations in which atoms continually move into and out of the FLL. In the FK model, the particles experience the substrate through a potential with a simple sinusoidal variation. This is to be compared with a substrate potential due to the summation of the Lennard-Jones potentials of the solid substrate atoms. Not only is the Lennard-Jones potential itself more complex than that of the sinusoidal FK potential, but the potential experienced by the liquid atoms further varies as they fluctuate along their $y$ and $z$ degrees of freedom. Additionally, the MD simulation incorporates finite temperature and heat fluxes that are absent in the FK model. Given these differences, the rather striking degree of similarity delivered by the FK model can be taken as evidence of the importance of geometry, particularly the ratio $r_{0}^{LL}/\lambda$, in governing the dynamics.

Finally, we note that unlike the coordinated kink motion described above, atomic motion into an unoccupied cell does not create the same extended atomic profile as seen in figure 12 and is not described by an FK kink. This difference is due to the asymmetry of the Lennard-Jones potential in the modified FK model (5.3). The Lennard-Jones potential produces the nonlinear wave profile of a kink because its nonlinear repulsive force is huge for liquid atoms brought close together. In contrast, it produces no such profile in the neighbourhood of a vacant cell, because it is only slightly attractive for liquid atoms separated across a vacant cell.

6. Conclusions

When sheared, liquid atoms move across a solid substrate in a motion known as slip. Slip has often been assumed to be diffusive, as liquid atoms move from one low-energy substrate site to another, independent of each other. However, there are conditions when slip takes place differently. Groups of atoms form into localized nonlinear waves that propagate at great speeds over the substrate – orders of magnitude faster than the slip velocity.

In this work, we focused on the atomic-level mechanisms by which liquid slip occurs, by systematically varying substrate lattice spacing $\lambda$ and characterizing the mechanisms of slip as a function of the parameter $r_{0}^{LL}/\lambda$. We have shown that there is a slip velocity minimum at $r_{0}^{LL}/\lambda =1$. Furthermore, under certain conditions, such as for $r_{0}^{LL}/\lambda <1$, slip occurs predominantly due to wave propagation. That is, the mass transport due to slip is conveyed predominantly by localized nonlinear waves, and not through surface diffusion. We have further shown that the observed waves are well described in their speed and profile by solutions to a modified Frenkel–Kontorova (FK) model and its continuum approximation, the sine-Gordon equation. Conversely, when $r_{0}^{LL}/\lambda >1$, liquid slip occurs predominantly due to mass propagation by individual atoms moving into unoccupied substrate sites, i.e. through surface diffusion. Slip under all circumstances is due to the sum of the contributions from these waves and the isolated motions of atoms into unoccupied substrate sites.

We have developed a novel conceptual framework to arrive at these conclusions. Using the doubly occupied cell as a numerical marker, we could identify individual localized propagating kinks in the noisy MD data. Fast 10 fs sampling was used to track kinks, revealing their well-defined profiles and fast speeds, which are comparable to predictions from the FK equation modified to use Lennard-Jones interactions between the liquid atoms.

The limitations of this work should be noted. The MD simulations utilize simple spherically symmetric Lennard-Jones atoms whose potential does not represent the complex potentials of even simple diatomic molecules. The substrate chosen for this research is also simple in structure and symmetric in its potential. Other crystal structures and real surfaces with defects present a complexity of substrate potential which is not modelled adequately here. The success of kinks to make their way through these more complex and realistic conditions remains to be investigated. Furthermore, as seen in figure 8, kinks are most prominent for $r_{0}^{LL}/\lambda <1$. Very low values of $r_{0}^{LL}/\lambda$ are problematic as the spacing, as measured by $\lambda$, of the substrate atoms can become so large as to allow for interpenetration of the liquid atoms. For neon over graphene, FCC gold, and diamond–cubic silicon substrates, $r_{0}^{LL}/\lambda$ is equal to 1.27, 1.08 and 0.81, respectively (Martini et al. 2006; Kannam et al. 2011; Celebi & Beskok 2018). Whether conditions that favour kinks will be easy to achieve in practice remains a question.

Solid-on-solid friction has been modelled using two opposing periodic lattices. As one lattice is forced over the other, strain deforms the lattices, leading to a misfit between the wavelengths of the two lattices that propagates as a kink (Braun et al. 2001; Woods et al. 2014). Sisan & Lichter (2014) considered a nanotube composed of hexagonal rings of carbon atoms containing a single file of TIP3P water molecules. With an applied pressure gradient, episodic pulses of flow occurred due to FK solitons that propagated along the length of the nanotube. These very different scenarios share with the present work two unequal length scales whose competition yields propagating kinks. These related researches with different materials and with various substrate potentials suggest that FK-type kinks may contribute to slip in settings more general than the cubic lattice structure studied here.

Perhaps there are practical applications of the new understanding of the slip presented here. Patterning the substrate to promote kinks may enhance slip or facilitate mixing at small scales. At a moving contact line, the inherent non-uniformity of relative atomic positions and speed of motion might be a meaningful missing contribution in treating the microscopic origins of contact line motion, microbubbles, and the nucleation of cavitation. The inhomogeneous distribution of momentum transfer at the small scales investigated here may also produce variations of heat transfer rate at the surface that could be exploited usefully. As a final speculation, we note that liquid–solid interfaces are home to commercially important catalysis. The kinetics of some reactions are found to confound typical descriptions in terms of surface diffusion. Furthermore, the catalytic dynamics on substrates that are densely populated by reactants and products is poorly understood. Perhaps the new understanding of atomic motion at the liquid–solid interface will provide a useful point of view to help to explain the complex cooperative events observed during catalysis on surfaces (Henß et al. 2019).