Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Jaffar
Hasnain
*^{a},
Swetlana
Jungblut
^{a},
Andreas
Tröster
^{b} and
Christoph
Dellago
^{a}
^{a}Faculty of Physics, University of Vienna, Austria. E-mail: jaffar.hasnain@univie.ac.at
^{b}Institute 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.

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.

(1) |

The substrate force, F_{sub}(r_{i}), acting on particle i is the negative gradient of the externally applied potential, U_{sub}(r_{i}) = −(U_{0}/9){3 + 2[cos(k_{1}r_{i}) + cos(k_{2}r_{i}) + cos(k_{3}r_{i})]}, where the k-vectors are chosen from the set k_{i}/‖k‖ ∈ { (0, 1)} with norm . This choice of the k-vectors produces a hexagonal arrangement of potential wells with a lattice constant a and lattice vectors . The depth of the potential minima was set to U_{0} = 27k_{B}T, where k_{B} 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 F_{max} = 24πk_{B}T/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 U_{Yuk}(r) = e^{−κr}/r, where the inverse screening length, κ = 6.25 μm^{−1}, determines the range of the Yukawa interaction. The coupling parameter is related to the effective charge on each colloid. Instead of specifying , we define the interaction strength as Γ = 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, F^{i}_{random}, with zero mean and variance 2k_{B}Tγ. The friction constant γ is related to the diffusion constant of a single particle by the Einstein relation, γ = D/k_{B}T. 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, F_{d} = (F_{d}, 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 × 10^{6} time steps for a multitude of values of Γ and F_{d} while keeping κ, a, and F_{max} constant. Before performing measurements, the systems were equilibrated for 10^{5} time steps. We used the reduced units γ = k_{B}T = 1 and all distances were rescaled by the lattice constant a and all forces by F_{max}.

In Fig. 2, this process is illustrated for a monolayer driven at two slightly different values of F_{d}. 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”.

Fig. 2 Motion of a monolayer with Γ/k_{B}T = 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, σ_{x}^{2} (third row), and the net force acting on the monolayer, F_{net}/F_{max} (bottom row), all as a function of reduced time, γt, for two values of F_{d}/F_{max}. 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, σ_{x}^{2}, 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^{−4}a^{2}, 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 F_{d}/F_{max} = 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 F_{d}/F_{max} = 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 F_{d}, yet the monolayer glides unhindered in the driving direction. In the following section, we will explain the origin of this “thermal sliding”.

(2) |

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

(3) |

(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 :

(4) |

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, u_{i}. We can therefore rewrite the position of each particle as r_{i} = u_{i} + R + R_{i}, where R_{i} 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, . For small u_{i}, 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 u_{i} is

(5) |

(6) |

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

(7) |

(8) |

(9) |

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, F^{eff}, is an (essentially) intensive quantity. The functional form of F^{eff}(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 F_{max}. 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 σ_{y}^{2}(R) when Γ/k_{B}T = 0.804, for which the average distance between a data point and the curve is 13% and the corresponding error in the estimation of F^{eff}(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).

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 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, t_{total}, that the monolayer takes to travel forwards by one lattice constant is an exponentially modified Gaussian distribution arising from the convolution of these distributions:

(10) |

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 F_{d} 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 F^{eff}(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 = [F_{d}^{2} − F^{eff}_{max}^{2}(Γ)]^{1/2}. The lines in Fig. 5 are plots of this simplification, and conform surprisingly well with the simulation results for F_{d} > F^{eff}_{max}(Γ), 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 F_{d} = F^{eff}_{max}, where velocities tend to be very small, then the obtained data would suggest the existence of a static friction located at F_{d} = F^{eff}_{max}(Γ). This apparent static friction obeys Amontons' law, in that the value of F^{eff}_{max}(Γ) 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 F_{max} ∝ U_{0}.^{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 F_{d} = F^{eff}_{max} 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}

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 Γ/k_{B}T = 1.038 as a function of F_{d} 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 F^{eff}_{max}/F_{max} = 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 Γ/k_{B}T = 0.4 produce both hopping wave and defect induced motion. We summarize our findings in Fig. 6, where we have plotted F^{eff}_{max}(Γ) for various values of aF_{max}/Γ (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.

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.

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

## Footnote |

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

This journal is © The Royal Society of Chemistry 2014 |