Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Frictional dynamics of stiff monolayers: from nucleation dynamics to thermal sliding

Jaffar Hasnain *a, Swetlana Jungblut a, Andreas Tröster b and Christoph Dellago a
aFaculty of Physics, University of Vienna, Austria. E-mail:
bInstitute of Theoretical Physics, Vienna University of Technology, Austria

Received 2nd April 2014 , Accepted 17th June 2014

First published on 20th June 2014

The inherently nonlinear dynamics of two surfaces as they are driven past each other, a phenomenon known as dry friction, has yet to be fully understood on an atomistic level. New experiments on colloidal monolayers forced over laser-generated substrates now offer the opportunity to investigate friction with single-particle resolution. Here, we use analytical theory and computer simulations to study the effect of thermal fluctuations on the stick-slip mechanism characteristic for the frictional response of a stiff colloidal monolayer on a commensurate substrate. By performing a harmonic expansion of the energy and employing elementary statistical mechanics, we map the motion of the monolayer onto a simple differential equation. Analytical expressions derived from our approach predict a transition from nucleation dynamics, where the monolayer moves in a sequence of activated hops over energy barriers, to “thermal sliding”, in which the effective substrate barrier opposing the motion of the monolayer disappears due to thermal fluctuations, leading to continuous, uninterrupted sliding motion. Furthermore, we find that the average velocity of the monolayer for large driving forces obeys a simple scaling behavior that is consistent with the existence of a static friction. For small forces, however, nucleation provides a mode of motion that leads to a small but non-vanishing mobility of the monolayer. Data obtained from simulations confirm this picture and agree quantitatively with our analytical formulae. The theory developed here holds under general conditions for sufficiently strong inter-particle repulsions and it yields specific predictions that can be tested in experiments.

1 Introduction

Although macroscopic laws of friction are centuries old and sufficiently accurate for a multitude of applications,1 and mesoscopic treatments have yielded great insights into this phenomenon,2,3 advances in nanotechnology require an understanding of friction on length and time scales in which atomistic details become important.4,5 Experimentally, several recent approaches, including atomic force microscopy setups,6–9 quartz-microbalance setups,10–13 and monolayers of charge-stabilized colloidal particles exposed to light-induced interference patterns,14 have made it possible to study friction on an atomistic level. In all of these cases, insights are gleaned from the dynamical response of small samples of particles that are driven across an external potential landscape. In particular, the remarkable accuracy and level of control achieved in experiments using two-dimensional crystals in laser fields14 have shed new light on the dynamics of collective excitations lying at the origin of friction. These experiments have been complemented by simulations and analytical studies of the Frenkel–Kontorova (FK) model.15–25 Despite its simplicity, this model captures the rich dynamical behavior and, more specifically, the soliton and anti-soliton structures largely determining the magnitude of friction observed in experiments.

Here, our aim is to investigate the effect of thermal fluctuations on the frictional response of stiff monolayers, i.e., monolayers with large elastic moduli, driven over a commensurate substrate. In a previous work,26 we examined the friction induced by a substrate with potential wells arranged in a regular structure identical to that of the unperturbed monolayer driven over it. In particular, we were interested in the role that the inter-particle interaction strength had on the sliding phases that the monolayer adopted as it was driven over such a commensurate substrate. We found that the inter-particle interaction strength regulates how many defects appear in the monolayer as it slides. When the interaction strength is so large that no defects appear in the system, the monolayer adopts a type of stick-slip motion which we call a “hopping wave”. This sliding mechanism consists of long periods of virtual motionlessness which are interrupted by the formation of small nuclei of particles that have escaped from their respective substrate wells. These particles subsequently initiate a cascade of particle hops which eventually encompass the entire system (Fig. 1). The “hopping wave” mechanism is of particular note, since it accounts for a sliding phase in which the monolayer remains structurally intact. It also appears to be a feature of systems of this type because the same mechanism was observed under different conditions in an earlier simulation study conducted by Reguzzoni et al.17 The purpose of our present work is to provide a general and quantitative understanding of the formation of these hopping waves.

image file: c4nr01790k-f1.tif
Fig. 1 Snapshot of a stiff monolayer driven by a constant force, Fd, that acts from top to bottom. Each of the particles is colored according to its substrate potential value. Red particles are at the top of their substrate potential barriers, purple particles are at the bottom of the wells, and green particles are somewhere along the walls.

2 Model

We simulated a two-dimensional array of overdamped particles exposed to an external substrate potential while being driven by a constant force. The particles repel each other via Yukawa interactions and receive random kicks due to the solvent surrounding them. The equation of motion for particle i, with position vector ri, is therefore
image file: c4nr01790k-t1.tif(1)

The substrate force, Fsub(ri), acting on particle i is the negative gradient of the externally applied potential, Usub(ri) = −(U0/9){3 + 2[cos(k1ri) + cos(k2ri) + cos(k3ri)]}, where the k-vectors are chosen from the set ki/‖k‖ ∈ {image file: c4nr01790k-t2.tifimage file: c4nr01790k-t3.tif (0, 1)} with norm image file: c4nr01790k-t4.tif. This choice of the k-vectors produces a hexagonal arrangement of potential wells with a lattice constant a and lattice vectors image file: c4nr01790k-t5.tif. The depth of the potential minima was set to U0 = 27kBT, where kB is the Boltzmann constant and T is the temperature, and the lattice constant of the substrate was set to a = 6 μm, so that the filling fraction is one, i.e., the number of substrate minima matches the number of particles in the system. As a result of the choice of these parameters, the maximum force that the substrate can exert on a particle in the x direction is Fmax = 24πkBT/a. The colloidal particles we study are charge-stabilized and therefore repel each other via Yukawa interactions. Two particles separated by a distance r from each other have a potential energy of UYuk(r) = [capital Gamma, Greek, tilde]eκr/r, where the inverse screening length, κ = 6.25 μm−1, determines the range of the Yukawa interaction. The coupling parameter [capital Gamma, Greek, tilde] is related to the effective charge on each colloid. Instead of specifying [capital Gamma, Greek, tilde], we define the interaction strength as Γ = [capital Gamma, Greek, tilde]eκa/a, which is the potential energy of two particles separated by one lattice constant. The choice of the parameter values used here resulted in the best agreement between simulation and experiment.14,26 For particles obeying overdamped Langevin dynamics, at each moment in time a particle experiences an uncorrelated Gaussian random force, Firandom, with zero mean and variance 2kB. The friction constant γ is related to the diffusion constant of a single particle by the Einstein relation, γ = D/kBT. We prepared a rectangular simulation box with periodic boundary conditions compatible with N = 6400 hexagonally arranged substrate minima. At the beginning of each simulation run, we placed a single particle in each well and applied a constant driving force, Fd = (Fd, 0), in the x direction. The algorithm presented in ref. 27 was employed to simulate the motion of the particles with a time discretization of δt = 10−4. Data were gathered from 100 runs of 2 × 106 time steps for a multitude of values of Γ and Fd while keeping κ, a, and Fmax constant. Before performing measurements, the systems were equilibrated for 105 time steps. We used the reduced units γ = kBT = 1 and all distances were rescaled by the lattice constant a and all forces by Fmax.

3 Buildup phase

The stick-slip mechanism mentioned above appears when the inter-particle interaction strength Γ is large enough, in comparison to the substrate potential depth, to preserve the hexagonal structure of the monolayer at all times. The monolayer, if it is to slide, has no other recourse than to form a distortion wave that travels through the system until each particle has moved forwards by one lattice constant. This is accomplished by forming a localized cluster of particles that escape their respective substrate wells. This cluster of hopping particles is the source of a sequence of particle hops that travels through the entire system. After such a “hopping wave” has run its course, the monolayer re-equilibrates and awaits the formation of yet another cluster (see ESI videos V1 and V2).

In Fig. 2, this process is illustrated for a monolayer driven at two slightly different values of Fd. In the top row of the figure, the average displacement of each particle from its position at t = 0 is plotted as a function of time. One can clearly see that the trajectories alternate between a “buildup phase”, in which the monolayer is in the process of forming a hopping wave, and a “slipping phase”, in which a hopping wave travels through the system until each particle has moved forwards by one lattice constant. Although the applied driving forces differ by only a fraction of a percent, the velocity of the monolayer changes almost by a factor 8. Moreover, the regularity with which the hopping waves appear changes drastically. The plateaus in the graph are colored in blue and correspond to data points belonging to the “buildup phase”.

image file: c4nr01790k-f2.tif
Fig. 2 Motion of a monolayer with Γ/kBT = 0.804. We have depicted the absolute displacement of the system, d (first row), the “periodic position of the center of mass”, R/a (second row), the variance in the x direction of the particle positions from the position of the center of mass, σx2 (third row), and the net force acting on the monolayer, Fnet/Fmax (bottom row), all as a function of reduced time, γt, for two values of Fd/Fmax. Data points colored in blue indicate that the monolayer is in the “buildup phase” whereas data points are colored red when a hopping wave is traveling through the system. Panels in the same column belong to the same trajectories. For animations of the data, see ESI videos V1 and V2.

We can take advantage of the periodicity of the substrate by considering the “periodic position of the center of mass” R (second row of Fig. 2), which is defined as the average displacement in the x direction of each particle from its nearest substrate potential minimum. During the buildup phase, the value of R is equal to the average displacement of the monolayer, d, modulo the lattice constant a. During the slipping phase, R does not have a physically meaningful value but the sharp peaks are nonetheless clear indicators of the existence of a hopping wave. The value R = 0.25a is of particular significance since it is the point of maximum resistance of the substrate. The plots of R vs. γt not only reveal the repetitive nature of the hopping wave mechanism, but also show that the monolayer on the left gets pinned by the substrate and takes a variable amount of time to evolve a hopping wave, whereas the buildup phase of the monolayer on the right consists of an essentially continuous drift towards the point of maximum resistance of the substrate, and shortly thereafter forms a hopping wave.

One can measure how accurately R describes the positions of the particles in the monolayer by calculating the variance in x direction of the particles' positions, σx2, from the periodic center of mass (third row of Fig. 2). As can be expected, during the buildup phase, the particles perform small oscillations about R whereas when a hopping wave is traveling through the system, large deviations are observed. We therefore consider the monolayer to be in the buildup phase if the variance is below a cutoff of 8.5 × 10−4a2, which evidently captures the plateaus in the first rows of the columns in Fig. 2.

An examination of the net force acting on the monolayer (bottom row of Fig. 2) confirms that the monolayer driven at Fd/Fmax = 0.980 experiences zero net force for significant periods of time, before a random fluctuation creates a hopping wave. The monolayer driven at a rate of Fd/Fmax = 0.984, on the other hand, is perpetually in motion. This is surprising since the substrate potential is clearly capable of applying a larger restoring force on each particle than Fd, yet the monolayer glides unhindered in the driving direction. In the following section, we will explain the origin of this “thermal sliding”.

4 Mapping onto the harmonic crystal

The trajectories in Fig. 2 differ primarily in the character and duration of their buildup phases. It also turns out that this part of the trajectories is amenable to analytical treatment. Consider the equation of motion of the periodic center of mass of the monolayer during the buildup phase,
image file: c4nr01790k-t6.tif(2)
where the stochastic force [F with combining macron]random(t) is the average random force acting on the particles and is a Gaussian uncorrelated noise with zero mean and variance 2kB/N. The Yukawa forces cancel due to Newton's third law and so we obtain a stochastic differential equation describing the motion of a Brownian diffuser in a potential defined by the driving force and the average substrate force acting on the monolayer. Since there is no diffusion in the direction perpendicular to the driving force (in the y direction), we merely need to consider the x component of the equation of motion of the system. The only unknown in the equation is the average substrate force acting on the monolayer in the x direction, image file: c4nr01790k-t7.tif, which depends on the distribution of the particles' positions ri and the instantaneous position of the center of mass, R.

In the following, we consider an ensemble of various realizations of configurations with a given R, as they appear in the system at different times along different trajectories. It then follows that the mean velocity of the monolayer during the buildup phase is

image file: c4nr01790k-t8.tif(3)
where the average 〈…〉R is performed over multiple buildup phases and Feff(R) is the effective force that acts on the center mass of the monolayer. To find an analytical expression for Feff(R), we propose the following two assumptions:

(1) During the buildup phase, the monolayer is in quasi-static thermodynamic equilibrium. It can be expected that this assumption holds because the monolayer creeps along the external potential landscape very slowly during the buildup phase.

(2) The total potential energy of the system can be approximated by a second order Taylor expansion. This approximation is expected to be valid for stiff crystals.

If the first condition is met, then, for a given position R, the probability distribution of the particle positions is the equilibrium distribution with the restriction image file: c4nr01790k-t9.tif:

image file: c4nr01790k-t10.tif(4)
where β = 1/kBT, δ is the Dirac delta function, and [r with combining right harpoon above (vector)] is a 2N-dimensional vector of x and y coordinates of each particle. Utot([r with combining right harpoon above (vector)]) denotes the total potential energy of the system as a function of all particle positions, and the normalization is given by the 2N-dimensional integral, image file: c4nr01790k-t11.tif. The effective substrate force, Feff(R), defined in eqn (3), is the corresponding ensemble average. The normalization and expectation values of such distributions typically defy analytical evaluation, so a further approximation is required.

The second assumption implies that the position of a particle deviates from the location of the center of mass, R = (R, 0), only by a small displacement, ui. We can therefore rewrite the position of each particle as ri = ui + R + Ri, where Ri is the position vector of the substrate minimum closest to the particle i. The center of mass of the monolayer is able to oscillate in the y direction but, by construction, the deviations in the x direction must cancel, image file: c4nr01790k-t12.tif. For small ui, the total potential of the system can then be approximated by a second order Taylor expansion. Although conceptually simple, the following calculation is rather cumbersome, thus, we have detailed every step of the procedure in the ESI Appendix S1 and mention only the essential results in the following. The work by Baumgartl et al.28 on this approach is recommended for further reading, and an introduction to the harmonic crystal can be found in ref. 29. The second order Taylor expansion of the total potential energy for small ui is

image file: c4nr01790k-t13.tif(5)
where the 2N × 2N interaction matrix [Doublestruck D](R) encodes how each particle interacts with the substrate and with all other particles and [u with combining right harpoon above (vector)]T is the transpose of the 2N-dimensional vector of the small displacements of the particle positions. Both in simulation and in this calculation, nearest neighbor cutoffs were used and, hence, the majority of the entries in [Doublestruck D](R) are zero. The function G(R) combines the linear term in the Taylor expansion of the substrate potential and the driving force Fd, and H(R) is the potential energy of the system if all ui are zero. Since the interaction matrix is symmetric, it can be brought to a diagonal form. Using the discrete Fourier transformation, the 2N × 2N matrix [Doublestruck D](R) can be resolved into N 2 × 2 matrices.28 The total energy of the system can therefore be rewritten as
image file: c4nr01790k-t14.tif(6)
where ū(q) and image file: c4nr01790k-t20.tif(q, R) are the discrete Fourier transforms of u and [Doublestruck D](R), respectively, and q is an element of the set of N reciprocal lattice vectors in the Brillouin zone of the substrate lattice. It turns out that the sum over all uix appearing in eqn (5) is image file: c4nr01790k-t15.tif, that is, the x component of the Fourier amplitude assigned to q = 0.

Recasting the total potential in this form allows us to evaluate Z(R), the expectation values of the variances of the particles' positions, and the average force that the substrate exerts on the monolayer for a given value of R. The expectation value of the force exerted by the substrate in the x direction is

image file: c4nr01790k-t16.tif(7)
where ∫d[u with combining right harpoon above (vector)] and ∫dū([q with combining right harpoon above (vector)]) are 2N-dimensional integrals over all independent degrees of freedom of the system. The only coordinate appearing in the delta function is ūx(0) and therefore the corresponding integral over this degree of freedom is unity. Since the linear term G(R) is only coupled to ūx(0), the result of the integration does not depend on G(R). The function H(R) can be pulled out of the integral and cancels with the corresponding term appearing in Z(R). What remains in eqn (7) is the product of 2N − 1 Gaussian integrals that can be evaluated individually for each degree of freedom. Having calculated the variances of ux, uy, and their cross correlation, we arrive at the result:
image file: c4nr01790k-t17.tif(8)
where σx2(R) and σy2(R) are the expectation values of the variances of the particles from their mean position in x and y directions. It turns out that, in the harmonic approximation of the hexagonal lattice, the cross correlation σxy vanishes which is compatible with our findings from simulation. These variances are directly related to the diagonal elements of the inverse of the dynamical matrix image file: c4nr01790k-t21.tif(q):
image file: c4nr01790k-t18.tif(9)
where μ ∈ {x, y}. The Kronecker delta δ term arises from the fact that the monolayer is free to oscillate in the y direction and ensures that an additional term is added to the sum when σy2(R) is calculated.

We find that the expectation values of the effective substrate force and the variances in the particles' displacements are independent of the applied driving force. Furthermore, the formulae for the variances in eqn (9) can be interpreted as the average value of a function and therefore the effective resistance due to the substrate, Feff, is an (essentially) intensive quantity. The functional form of Feff(R) is similar to the original external potential except that it is exponentially reduced in terms of the variances, σμ2(R), and thus the temperature T. For a sinusoidal substrate, each of the variances increases monotonically as the monolayer travels along the barrier. We included an instruction file, code for a C program, a python script, and a pair of sample files in the ESI Appendices S3–S7, with which these formulae can be evaluated.

In Fig. 3, we compare the curves generated by eqn (8) and (9) with data obtained from simulation during the buildup phase. We used 4 monolayers of varying interaction strengths Γ, and applied various driving forces close to, but less than Fmax. Data from simulations were compiled by calculating the variances of the particle positions and the net substrate force acting on the center of mass, all as a function of the center of mass R, during the buildup phase. The mean values of these quantities were then plotted along with the theoretical predictions in Fig. 3. A more detailed exposition of how the data were gathered and analyzed can be found in the first section of the ESI Appendix S2. As predicted by our calculations, there is a range of R/a in which the driving force does not influence the expectation values of the monolayer, which explains the collapse of the data points of the same color within the indicated region. Furthermore, the lines, which are our theoretical predictions, conform very well with the simulation data. We consider the main source of error to be the truncation of the Taylor expansion at the second order because the theoretical curves become more accurate as the interaction strength is increased. The largest error is incurred in the estimation of σy2(R) when Γ/kBT = 0.804, for which the average distance between a data point and the curve is 13% and the corresponding error in the estimation of Feff(R) is 6.6%. Although the region of space in which our formulae are valid may seem small, the monolayers reside in this region for the vast majority of the time. The harmonic approximation diverges shortly after R = 0.25a because the curvature of the external potential landscape changes sign and therefore ceases to be a pinning potential. This change of curvature explains why, in simulation, hopping waves form almost immediately after the monolayer reaches that point (see right panels of Fig. 2).

image file: c4nr01790k-f3.tif
Fig. 3 Comparison of theoretical prediction (lines) with data gathered from simulation runs (symbols). Symbols and lines of the same color correspond to the same value of Γ, whereas different symbols of the same color correspond to different driving forces. The dashed black line in the bottom panel of the graph is Feff(R) for T = 0 or, equivalently, Γ = ∞. Vertical black lines limit the region in which the presented theory is applicable.

5 Dynamical phases

In the previous section, we mapped the buildup phase of the motion of the monolayer onto the motion of a single overdamped Brownian diffuser subject to an effective substrate force, Feff(R) (eqn (2) and (3)). This effective substrate force is necessarily weaker than the substrate force acting on each particle due to the thermal motion of the particles in the monolayer. The formula for Feff(R) can be used to find the value of the largest restoring force due to the effective substrate, Feffmax, which delimits two different dynamical regimes. In the first regime, when Fd > Feffmax, the monolayer slides unhindered over an effective substrate. Hence we call this motion “thermal sliding”. In the second regime, when the driving force is lower than the effective barrier, Fd < Feffmax, the monolayer becomes pinned and its center of mass fluctuates about an equilibrium position which can be inferred from Fig. 3. In this latter case, a critical number of particles need to be kicked out of their respective potential wells in order to initiate a hopping wave. The time taken to form such a critical cluster needs to be treated within the framework of the classical nucleation theory.17

The entire trajectory of a monolayer can therefore be resolved into three phases: the buildup phase, the nucleation phase (where applicable), and the hopping wave phase. We expect that the distributions of both the time taken to complete the buildup phase and the time for a hopping wave to travel through the system are Gaussians. Since the convolution of two Gaussian distributions is also a Gaussian, we define [t with combining tilde] to be the mean time that the monolayer takes to complete the buildup phase plus the mean time the hopping wave takes to travel through the system, and the quantity ζ2 is the sum of the variances of the aforementioned times. The time taken for nucleation to occur, on the other hand, obeys an exponential distribution with a characteristic time τ. The distribution of the total time, ttotal, that the monolayer takes to travel forwards by one lattice constant is an exponentially modified Gaussian distribution arising from the convolution of these distributions:

image file: c4nr01790k-t19.tif(10)
where erfc(x) is the complementary error function and Π is the product of the norms of the distributions. In the top panel of Fig. 4, we have plotted the distribution of the waiting times between two successive hopping waves, P(ttotal), for a monolayer driven with a force less than, equal to, and greater than Feffmax. The data were gathered by finding the time between two successive peaks in the net force acting on the monolayer (see bottom row of Fig. 2). The solid lines are fits of eqn (10) to the data sets, and evidently the tail of the distribution (determined by the value of τ) disappears when Fd becomes larger than Feffmax. In the bottom of Fig. 4, we plotted the mean nucleation time, τ, gathered from fits of eqn (10), for different monolayers as a function of the driving force. Due to the exponential dependence of τ on the driving force, and the fact that the nucleation barrier is expected to vanish if Fd > Feffmax, one can observe a change of at least 1 order of magnitude in the nucleation time within a window of 0.4% of Feffmax. The shape of the curves in the bottom of Fig. 4 indicates that the dynamical transition from nucleation to thermal sliding is continuous.

image file: c4nr01790k-f4.tif
Fig. 4 Top: distribution of waiting times between hopping waves for a monolayer with Γ/kBT = 0.804 for driving forces below, equal to, and above the effective barrier height, Feffmax/Fmax = 0.982. The red and black curves have been shifted horizontally by 15γt and 40γt, respectively, for clarity. Bottom: nucleation time, γτ, extracted from the distribution of waiting times for multiple monolayers driven at different rates.

In Fig. 5, we have plotted the mean velocity of a number of monolayers driven at rates both above and below their respective effective force barriers. All of the considerations herein lead to the expectation that there are two different scaling regimes of the mean velocity of the monolayer, γν, with respect to the driving force. In the nucleation regime, also referred to as the creep regime, the value of the mean nucleation time, τ, has the largest influence on the velocity of the monolayer and scales exponentially in the height of the free energy barrier posed by the substrate. The free energy barrier, in turn, is influenced by the driving force and thus the velocity is expected to decay exponentially as Fd goes to zero. In the thermal sliding regime, the mean velocity of the monolayer depends primarily on the time the center of mass takes to diffuse along the potential landscape associated with Feff(R). If one ignores the random force in eqn (2) and assumes that the variance in the particle positions remains constant, then the first order differential equation for the motion of the monolayer can be solved quite easily. The mean velocity of the monolayer for the thermal sliding regime is then given by γv = [Fd2Feffmax2(Γ)]1/2. The lines in Fig. 5 are plots of this simplification, and conform surprisingly well with the simulation results for Fd > Feffmax(Γ), especially with respect to the scaling of γν. As a corollary to this consideration, if, for some reason one were unable to measure velocities in the regime close to Fd = Feffmax, where velocities tend to be very small, then the obtained data would suggest the existence of a static friction located at Fd = Feffmax(Γ). This apparent static friction obeys Amontons' law, in that the value of Feffmax(Γ) is independent of the contact area, which in this model is the particle number N, and is, to the first order, proportional to the applied load, which is FmaxU0.1,30 Nonetheless, we have shown that the atomic details of such monolayers induce a dramatic change in the scaling of the mean velocity close to Fd = Feffmax and that it decays exponentially for small driving forces. This finding may have some bearing on the discussion of the molecular origin of static friction found in the literature.17,21,30–32

image file: c4nr01790k-f5.tif
Fig. 5 Velocity profiles of the simulated monolayers as a function of the reduced driving force Fd/Fmax. The lines are plots of the function γv = [Fd2Feffmax2(Γ)]1/2 that terminate at Fd = Feffmax(Γ) and roughly reproduce the velocity of the monolayers when Fd > Feffmax(Γ). Inset: mean velocity of a monolayer with Γ/kBT = 1.038 for different system sizes N.

According to nucleation theory, the mean time of observing a nucleation event is inversely proportional to the product of the nucleation rate and the volume of the system τ = [JV]−1. Therefore, when a collective mechanism is responsible for the sliding motion, one can expect that an increase of the system size will result in an increase of the monolayer mobility, contrary to the traditional experiences with friction. The mean velocity is expected to converge, however, for sufficiently large volumes or high nucleation rates, due to the appearance of multiple, simultaneous, hopping waves. The inset in Fig. 5 is a plot of the mean velocity of a monolayer with Γ/kBT = 1.038 as a function of Fd for different system sizes N, and, as one can clearly see, the mobility of the monolayer increases with N in the nucleation regime, which is expected to end at Feffmax/Fmax = 0.9855, and indeed, shortly thereafter, the velocities converge. A more detailed examination of the behavior of the system as a function of N can be found in the second section of the ESI Appendix S2.

We attempted to find a suitable criterion to determine under which conditions the hopping wave mechanism gives way to defect-driven motion, but we have been, so far, unsuccessful. We do expect, however, that the ratio of the substrate potential depth (which favors the formation of defects) and the interaction strength (which penalizes defects) is to the first order the main factor in determining whether the system moves through the formation of hopping waves or through correlated defects. Part of the difficulty in finding a cutoff between these two mechanisms stems from the fact that this is a continuous transition, as illustrated in Video S3, which shows that monolayers with Γ/kBT = 0.4 produce both hopping wave and defect induced motion. We summarize our findings in Fig. 6, where we have plotted Feffmax(Γ) for various values of aFmax/Γ (red dashed line) to delineate the nucleation regime from the thermal sliding regime. The symbols indicate interaction strengths and driving forces that we simulated. Empty boxes represent runs in which the net force acting on the monolayer was zero for significant periods of time, as is necessary for nucleation to occur. Filled circles, on the other hand, denote simulations in which the net force was positive virtually all the time, as expected for thermal sliding.

image file: c4nr01790k-f6.tif
Fig. 6 Proposed dynamical phase diagram of overdamped monolayers driven by a constant force. The dashed red line is the theoretical demarcation between thermal sliding and nucleation-induced motion, given by Feffmax(Γ). Empty boxes denote parameter values for which nucleation of hopping waves was found and filled circles correspond to simulation runs in which thermal sliding occurred. The green gradient indicates the region in which the hopping wave mechanism gives way to defect-induced motion.

6 Conclusions and outlook

In this work, we found that stiff overdamped monolayers, driven over commensurate substrates, adopt one of two different dynamic phases and we developed a quantitative theory to describe them. In particular, the velocity of the monolayer can be entirely characterized by the time the center of mass creeps along the effective substrate potential, the time a hopping wave travels through the system, and, where applicable, the mean nucleation time. While the results we presented in this work can be used to determine the first of the aforementioned times, analytical expressions for the hopping wave velocity and the nucleation time need to be developed. With these three quantities, the velocity profiles of the monolayers driven over a commensurate substrate can be reconstructed analytically. The driving force determines which of these times has the largest influence on the velocity of the system and a crossover occurs at Fd = Feffmax(Γ, Fmax, a, β), which has been computed analytically. For large driving forces, the mean velocity of this model scales as if there existed a static friction and, for small forces, we indicated that the mobility has an atypical dependence on the contact area, but a detailed analysis of this finite size effect is still pending.

The theory presented herein remains entirely unchanged if a different radially symmetric interaction potential is employed, provided that the particles in the system always repel each other strongly. We can therefore predict that a density dependent dynamical transition occurs for non-monotonic potentials even for large interaction strengths. The application of the approximation to substrates with different geometries, on the other hand, is much more difficult. It stands to reason, that if the geometry of the substrate is different from that of the monolayer, then the theory is only applicable if, during the buildup phase, each substrate minimum is occupied by exactly one particle. The work conducted by McDermott et al.22 and by Reguzzoni et al.17 are two examples of systems in which the shapes of the substrate and monolayer are very different, but the theory may remain applicable. Furthermore, although the dynamic phases of the more complex underdamped Langevin dynamics would introduce an additional parameter to the system, it ought to be manageable under the right conditions, particularly in the “onset of sliding” simulations performed in ref. 20.


This research was funded by the FWF project numbers P24681-N20, V 305-N27, and P22087-N16. All simulations were performed on the Vienna Scientific Cluster (VSC). The authors gratefully acknowledge Philipp Geiger, Georg Menzl, and Michael Grünwald for useful discussions.


  1. G. Amontons, Mem. Acad. Roy. Sci., 1699, 206–222 Search PubMed.
  2. N. Okamoto and M. Nakazawa, Int. J. Numer. Methods Eng., 1979, 14, 337 CrossRef.
  3. R. S. Sayles, Tribol. Int., 1996, 29, 639 CrossRef.
  4. M. Urbakh and R. Meyer, Nat. Mater., 2010, 9, 8 CrossRef CAS PubMed.
  5. Y. Mo, K. T. Turner and I. Szlufarska, Nature, 2009, 457, 1116 CrossRef CAS PubMed.
  6. B. Saha, E. Liu and S. B. Tor, Nanotribological phenomena, principles and mechanisms for MEMS, Springer-Verlag Berlin Heidelberg, 2013 Search PubMed.
  7. J. A. Ruan and B. Bhushan, J. Tribol., 1994, 116, 378 CrossRef CAS.
  8. M. Langer, M. Kisiel, R. Pawlak, F. Pellegrini, G. E. Santoro, R. Buzio, A. Gerbi, G. Balakrishnan, A. Baratoff, E. Tosatti and E. Meyer, Nat. Mater., 2014, 13, 173 CrossRef CAS PubMed.
  9. R. Pérez, Nat. Mater., 2014, 13, 118 CrossRef PubMed.
  10. J. Krim, D. H. Solina and R. Chiarello, Phys. Rev. Lett., 1991, 66, 181 CrossRef CAS.
  11. T. Coffey and J. Krim, Phys. Rev. Lett., 2005, 95, 076101 CrossRef CAS.
  12. L. Bruschi, A. Carlin and G. Mistura, Phys. Rev. Lett., 2002, 88, 046105 CrossRef CAS.
  13. J. Jupille, J. J. Ehrhardt, D. Fargues and A. Cassuto, Faraday Discuss. Chem. Soc., 1990, 89, 323 RSC.
  14. T. Bohlein, J. Mikhael and C. Bechinger, Nat. Mater., 2012, 11, 126 CrossRef CAS PubMed.
  15. Y. Frenkel and T. Kontorova, Phys. Z. Sowjetunion, 1938, 13, 1 CAS.
  16. L. Prandtl, Z. Angew. Math. Mech., 1928, 8, 85 CrossRef.
  17. M. Reguzzoni, M. Ferrario, S. Zapperi and M. C. Righi, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 1311 CrossRef CAS PubMed.
  18. O. M. Braun, A. R. Bishop and J. Röder, Phys. Rev. Lett., 1997, 79, 3692 CrossRef CAS.
  19. J. Tekić, O. M. Braun and B. Hu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 026104 CrossRef.
  20. O. M. Braun, M. V. Paliy, J. Röder and A. R. Bishop, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 036129 CrossRef CAS.
  21. A. Vanossi, N. Manini and E. Tosatti, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 16426 CrossRef PubMed.
  22. D. McDermott, J. Amelang, C. J. Olson Reichhardt and C. Reichhardt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062301 CrossRef CAS.
  23. C. Reichhardt and C. J. Olson Reichhardt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051401 CrossRef CAS.
  24. O. M. Braun and Y. S. Kivshar, The Frenkel–Kontorova model: concepts, methods, and applications, Springer-Verlag Berlin Heidelberg, 2004 Search PubMed.
  25. A. Vanossi and O. M. Braun, J. Phys.: Condens. Matter, 2007, 19, 305017 CrossRef.
  26. J. Hasnain, S. Jungblut and C. Dellago, Soft Matter, 2013, 9, 5867 RSC.
  27. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, 1989 Search PubMed.
  28. H. H. von Grünberg and J. Baumgartl, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 051406 CrossRef PubMed.
  29. N. W. Ashcroft and N. D. Mermin, Solid state physics, Brooks/Cole, Cengage Learning, 1976 Search PubMed.
  30. M. H. Müser, L. Wenning and M. O. Robbins, Phys. Rev. Lett., 2001, 86, 1295 CrossRef.
  31. M. H. Müser, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 1257 CrossRef PubMed.
  32. G. He, M. H. Müser and M. O. Robbins, Science, 1999, 284, 1650 CrossRef CAS.


Electronic supplementary information (ESI) available. See DOI: 10.1039/c4nr01790k

This journal is © The Royal Society of Chemistry 2014