{'fraction_non_alphanumeric': 0.03250263784795118, 'fraction_numerical': 0.011758819114335751, 'mean_word_length': 5.02122085407388, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Friction and Nonlinear Dynamics

The nonlinear dynamics associated with sliding friction forms a broad interdisciplinary research field that involves complex dynamical processes and patterns covering a broad range of time and length scales. Progress in experimental techniques and computational resources has stimulated the development of more refined and accurate mathematical and numerical models, capable of capturing many of the essentially nonlinear phenomena involved in friction.


Introduction
Friction phenomena take place across a broad range of time and length scales, from microscopic atomistic processes, as in the gliding motion of a nanocluster or a nanomotor [1], up to extremely macroscopic instances, as in fault dynamics and earthquake events [2]. Due to the ubiquitous nature of mechanical dissipative processes and the enormous practical relevance, friction has been investigated over the centuries. While the empirical laws of macroscopic friction are well known [3], the fundamental understanding of the tribological phenomena at the microscopic scales is still lacking from many points of view. The basic difficulty of friction is intrinsic, involving the dissipative dynamics of large systems, often across ill-characterized interfaces, and generally violent and nonlinear. The severity of the task is also related to the experimental difficulty to probe systems with many degrees of freedom under a forced spatial confinement, that leaves very limited access to probing the buried sliding interface. Thanks to remarkable developments in nanotechnology, new inroads are being pursued and new discoveries are being made. At the nanometer scale, state-of-the-art ultra-high-vacuum systems and local probe studies show a dynamical behavior which is often significantly different, not just quantitatively but even qualitatively, from the ones observed in macroscopic tribology. Bridging the gap among the different length scales in tribological systems still remains an open challenge. The phenomenological descriptions that apply to macroscopic friction cannot yet be derived from the fundamental atomic principles and the interplay of processes occurring at the molecular level. Nanofriction is in somewhat better shape. Together with the current experimental possibility to perform well-defined measurements on well-characterized materials at the fundamental microscopic level of investigation of the sliding contacts, advances in the computer modeling of interatomic interactions in materials science and complex systems encompass molecular-dynamics (MD) simulations of medium to large scale for the exploration of the tribo-dynamics with atomic resolution [4,5]. Despite the benefits brought about by numerical simulations of realistic 3D sliding systems, the resulting proliferation of detailed complex data, and the requirement of always-growing computational efforts have stimulated, in parallel, the concurrent search for simpler modeling schemes, such as, e.g., generalized Prandtl-Tomlinson (PT), Frenkel-Kontorova (FK), for nanofriction, and of Burridge-Knopoff and earthquake-like models, for mesoscale and macroscale friction, suitable to describe the essence of the physics involved in highly nonlinear and non-equilibrium tribological phenomena in a more immediate fashion.
Here we discuss current progress and open problems in the simulation and modeling of tribology at the microscopic scale, and its connection to the macroscale. Neither the PT model, described in detail in several surveys [5,6] with several applications to concrete tip-based physical systems, nor the phenomenological approach based on the rate-and-state models [7] will be considered here. With a view to emphasize the role of nonlinearity the present topical review will restrict to the following theoretical approaches to sliding friction. Section 2 revises the simple case of near-equilibrium linear friction in classical mechanics. Section 3 focuses on nonlinearity in crystal sliding in the framework of the FK model and its generalizations. Atomistic models and MD nanofriction simulations are presented in Sec. 4. Mesoscopic multicontact earthquakelike models are finally examined in Sec. 5.

Linear Friction and Dissipation
Statistical mechanics accounts for the intimate mechanism of friction: a system at equilibrium has its kinetic energy uniformly distributed among all its degrees of freedom. A sliding macroscopic object clearly is not at equilibrium: one of its degrees of freedom (the center-of-mass motion) has far more kinetic energy than any other. The tendency of the system toward equilibrium will lead to the transfer of energy from that degree of freedom to all other ones: as a result the macroscopic object will slow down and its energy will be transferred to the disordered motion of the other degrees of freedom, resulting in warming up. This is all sliding friction really is: the tendency of systems toward equilibrium energy equipartitioning among many interacting degrees of freedom.
Thus, in the course of friction under an applied external force, energy is reversed into the system in the form of frictional heat. The frictional heat is generally dissipated by some form of heat bath, such as that provided by a thermostat at temperature T . In a frictional steady state, caused for example by submitting a slider to an external force F , the slider dissipates energy to the bath, and therefore does not accelerate indefinitely -it reaches instead a steady state characterized by an average drift velocity v . When both v and F are infinitesimal, the relationship between the two quantities is linear, In this so-called "viscous friction" the proportionality constant γ is the linear friction coefficient. It is known from classical statistical mechanics, for example of Brownian motion as described by the Langevin equation, that for linear friction systems which obey Eq. (1), the Einstein relationship is generally valid, connecting the friction coefficient γ, which measures dissipation, to the diffusion coefficient D, which measures fluctuations. This expresses the fluctuationdissipation theorem of linear, viscous friction.

The Fluctuation-Dissipation Theorem
To simulate the classical motion of a macroscopic object moving in contact with an equilibrium bath such as the molecules of a gas or a liquid, or the phonons of a solid, the standard implementation [8,9] requires adapting Newton's equations of motion with the addition of a damping force f damp plus a random force f rand (t). The damping force represents the transfer of energy from the macroscopic object of mass m to the heat bath, i.e. dissipation: This formula, equivalent to Eq. (1), assumes that the deviation from equilibrium is small, so that linear response holds: the restoring Stokes force is linear in the perturbing velocity, and acts opposite to it to restore the ˙ r = 0 equilibrium regime. This linear dependence is purely the lowest order term in a Taylor expansion: there is no reason to expect the linear relation (3) to extend to large velocity, and indeed e.g. the drag friction of speeding objects in gases is well known to follow Rayleigh's quadratic dependence on speed |˙ r|. The random-force term represents statistically the "kicks" that the objects experiences due to its interaction with the thermal bath. In the frame of reference of the thermal bath f rand (t) = 0 of course. The random term is the result of many very frequent collisions events, resulting in random forces uncorrelated with themselves except over very short time spans. More precisely, we assume there is some maximum time τ beyond which any correlation vanish: In addition, the assumption of thermal equilibrium ensures us that the bath is in a steady state, so that f rand α (t)f rand α ′ (t + δt) is independent of t, and depends on δt only: the statistical properties of the random force are constant in time. In most practical situations, one can safely ignore the dynamics over a time scale of the order of τ or shorter. We are interested instead in the integral effect of f rand (t) over some time period t that is long compared to τ . We can break up that integral into many pieces, each covering a duration τ : f rand α (t ′ )dt ′ +. . . (5) This integral is then a sum of many independent random terms, each drawn from the same distribution whose only relevant property is that it has zero mean value. As a result of the central-limit theorem, the total integral obeys a Gaussian distribution with null mean, and whose standard deviation scales with the number of terms in the sum, i.e. t 1/2 . By taking the equation of motion in the absence of any external driving, and integrating it in time we obtaiṅ The first term at the right-hand side of Eq. (7) becomes negligible for a time t ≫ 1/γ, long enough for the object to equilibrate with the thermostat and lose memory of its initial condition. In this large-t limit, by taking the square module of Eq. (7) and executing the ensemble average, we have The term at the left side multiplied by 1 2 m yields the average kinetic energy, which by standard equipartition needs to equal 3k B T /2. By rearranging the exponentials on the right-hand side and substituting t 1 = t − t ′′ and t 2 = t ′′ − t ′ , we obtain: where we have used the steadiness of the stochastic process discussed after Eq. (4). Assuming, as is commonly the case, that the autocorrelation time τ of the random term is short compared to γ −1 , the integrand of Eq. (9) has e −γt 2 ≃ 1 in all region |t 2 | ≤ τ of delays where the factor f rand (0) · f rand (t 2 ) is significantly different from 0. This observation further simplifies Eq. (9) to the fluctuation-dissipation relation This expression draws an explicit link between the autocorrelation amplitude of the fluctuations and the product of the dissipation coefficient and thermostat temperature. Note that the integral in Eq. (10) depends on both the amplitude of fluctuations f rand and the time over which they remain self-correlated. The effect on the mesoscale object increases if the random force is larger and/or if the time interval over which f rand pushes in the same direction before changing is longer. The relation (10) can be equally well satisfied by weaker random forces acting for longer correlation times or stronger forces with shorter τ , which thus lead to the same statistical effects. As τ is the shortest time scale around, for all practical purpose one can satisfy Eq. (10) assuming a sort of τ → 0 limit: providing a simple recipe for computer simulations. For simulations of models such as the PT or FK ones, the phenomenological degrees of freedom are often coupled to a Langevin thermostat of this kind, implying that each degree of freedom is actually coupled to a vast number of other bath degrees of freedom. Even in MD simulations of atomic-scale friction, Langevin thermostats are applied to all or a part of the atoms involved [9][10][11]. Of course, this approach is not rigorous, since the relevant particles colliding with each given simulated atom are already all included in the conservative and deterministic forces explicitly accounted for by the "force field". The Langevin approach is quite accurate to describe small perturbations away from equilibrium, but it may fail quite badly in the strongly out-of-equilibrium nonlinear phenomena which are the target of the present paper.

Linear versus Nonlinear Friction
In the rest of this review we will deal with nonlinear frictional phenomena, which deviate violently from linearity and near-equilibrium, and where therefore Eqs. (1) and (2) do not generally apply. As it has, surprisingly, only been realized in the last few decades, even arbitrarily violently non-equilibrium and nonlinear driven phenomena adhere to an extension of the fluctuation-dissipation theorem. That is the Jarzynski (or Jarzynski-Crooks) relation [12,13], whose simplest form can be briefly summarized as follows. Suppose starting from a system in state A at temperature T , and apply an external force of arbitrary form and strength causing it to evolve, for example to slide, to another state B; assume for simplicity B to be also a state of equilibrium. Call W AB the work done by the external force, and call ∆F AB = F B − F A the difference of equilibrium free energy between the states B and A. Clearly, W AB > ∆F AB must be valid, because some work will be wasted in going from A to B, unless that was done infinitely slowly (adiabatically).
Suppose now repeating the forced motion A→B many times. Each time, W AB will be different. The Jarzynski equality states that It can be shown that in near-equilibrium conditions, Eq. (12) is completely equivalent to the fluctuation-dissipation theorem. The beauty of it is however that Eq. (12) is totally general.
One particular case is useful in order to underline its far-reaching power. Suppose to take B=A, that is a final state identical to the starting one. In that case This equation appears at first sight impossible to satisfy, because surely all W AA > 0: all forced motion must cost work. The answer is that the probability P (W AA ) is indeed a distribution centered around a positive W AA > 0, but with a nonzero tail extending to W AA < 0. This tail represents rare events where work is gained rather than spent -we can think of them as a sort of "free lunches". Jarzynski's theorem requires that "free lunches" must occur precisely in such a measure to satisfy Eq. (13). However it is easy to convince ourselves that they will be frequent only in microscopic systems, where P (W AA ) is broad. The larger the system involved, the narrower P (W AA ) will be, the rarer and rarer the occurrence of "free lunches". In a macroscopic friction experiment, the occurrence of a "free lunch" will be virtually impossible.

The Frenkel-Kontorova Model
In nanoscale tribology, extensive attention has focused on the time-honored PT model, which describes a point-like tip sliding over a space-periodic crystalline surface in a minimal fashion. We shall omit this model from the present review, since it is covered in great detail elsewhere [5,6,14]. We concentrate instead on its natural extension, the one-dimensional FK model [15], which provides a prototypical description of the mutual sliding of two perfect, extended crystalline surfaces. First studied analytically in Ref. [16] and later introduced independently to address the dynamics of dislocations in crystals [17][18][19], subsequently this model became the paradigm describing the structure and dynamics of adsorbed monolayers in the context of surface physics. The standard FK model consists of a 1D chain of N classical particles ("atoms"), interacting via harmonic forces and moving in a sinusoidal potential, as sketched in Fig. 1

. The Hamiltonian is
In Eq. (14), the p 2 i /(2m) term represents the kinetic energy of the particles, and the next term describes the harmonic interaction, with elastic constant K, of nearestneighboring atoms at equilibrium distance a c . The final cosine term describes the "substrate corrugation", i.e. the periodic potential of amplitude U 0 and period a b , as experienced by all particles alike. To probe static friction, all atoms are driven by an external force F , which is increased adiabatically until sliding starts.
The continuum limit of the FK model, appropriate for large K, is the exactly integrable sine-Gordon (SG) equation, and this mapping contributed to the great success of the FK model. The solutions of the SG model include nonlinear topological solitons (known as "kinks" and "antikinks"), plus dynamical solitons ("breathers"), beside linear vibration waves (phonons). In the FK model, the sliding processes are entirely governed by its topological excitations, the kinks. Let us consider the simplest "commensurate" case, where before sliding the chain is in a trivial ground state (GS), when N atoms fit one in each of the M minima of the substrate potential, so that the coverage (i.e. the relative atomic concentration) θ = N/M = a b /a c equals 1. In this case, the addition (or subtraction) of a single atom results in configurations of the chain characterized by one kink (or antikink) excitation. Still at zero applied force, in order to reach a local minimum of the total potential energy in Eq. (14), the kink expands in space over a finite length, so that the resulting relaxed chain configuration consists in a local compression (or expansion, for an antikink). Upon application of a force, it is far easier to move along the chain for kinks than for atoms, since the activation energy ε PN for a kink displacement [known as the Peierls-Nabarro (PN) barrier] is systematically smaller, and often much smaller [20], than the amplitude U 0 of the energy barrier that single atoms experience in the substrate corrugation.
The motion of kinks (antikinks), i.e. the displacement of the extra atoms (vacancies) represents the mechanism for mass transport along the chain. These displacements are responsible for the mobility, diffusivity, and conductivity within this model. Generally therefore a larger concentration of kinks is associated to a larger the overall mobility [22]. For the simple commensurate GS (e.g., θ = 1), which contains neither kinks nor antikinks, the onset of sliding motion under a driving force occurs via the creation of a kink-antikink pair, e.g. induced by a thermal fluctuation, see Fig. 2.
If the FK chain is of finite length, kinks/antikinks are usually created at one chain free end, then they advance along the chain, eventually disappearing at the opposite end [23]. Every kink running from one end of the chain to the other produces the advancement of the entire chain by one lattice spacing a b . For a finite film confined between two surfaces, or for an island deposited on a surface, the general expectation is that sliding initiates likewise with the formation and entrance of a kink, or antikink at the boundary [24]. In this two-dimensional (2D) case, and more generally in D-dimensional systems, the zero-dimensional kinks of the FK model are replaced by (D−1) dimensional misfit dislocations or domain walls, whose qualitative physics and role is essentially the same.
Incommensurability between the periods a b and a c plays an important role in the FK model. Assume, in the limit of an infinite chain length, the ratio θ = a b /a c of the substrate period a b to the average spacing a c of the chain to be irrational. The GS of the resulting incommensurate FK model is characterized by a sort of "staircase" deformation, with a regular sequence of regions where the chain is compressed (or expanded) to match the periodic potential, separated by kinks (or antikinks), where, at regular intervals, the misfit stress is released through a localized expansion (compression). The incommensurate FK model exhibits, under fairly general conditions on θ [25], a critical elastic constant K = K c , such that if K > K c the chain can slide freely on the substrate at no energy cost, i.e. the static friction F s drops to zero (and the low-velocity kinetic friction becomes extremely small), while, remarkably, this is no longer true when K < K c . In the early 1980s a rigorous mathematical theory of this phenomenon called "the transition by breaking of analyticity", now widely known as the Aubry transition [26][27][28][29][30][31], was developed. A simple explanation of free sliding in the unpinned F s = 0 state is the following. For every atom climbing up toward a corrugation potential maximum, there always is another atom moving down, with an exact energy balance of these processes. Quite generally, incommensurability guarantees that the total energy (we are at T = 0) is independent of the relative position of the chain and the periodic lattice. However, in order for the chain to slide with continuity between two consecutive positions, it is necessary that particles should be able to occupy a maximum of the potential, the worst possible position. At the Aubry transition, however, realized by a relative increase of the periodic potential magnitude, or equivalently by a softening of the chain stiffness, the probability for a particle to occupy that position drops from a finite value to exactly zero. The nature of this transition, which is structural but without any other static order parameter (besides energy, of course), is dynamical, similar in that to a glass transition: simply, a part of phase space becomes unavailable, in this case by sliding. The chain is unpinned and mobile as long as, in its GS, atoms may occupy with a finite probability all positions, including those arbitrarily close to the maxima of the substrate potential, but is immobilized when that possibility ceases. The critical chain stiffness K = K c marks the crossing of the Aubry transition, where the chain turns from the free sliding state to the locked ("pinned") state with a nonzero static friction F s [20]. The value K c is in turn a discontinuous function on the length ratio θ characterizing the model. The minimum value K c ≃ 1.0291926 [in units of 2U 0 (π/a b ) 2 ] is achieved for the golden-mean ratio θ = (1 + √ 5)/2 [15]. The stiff-spring chain with K > K c can explore adiabatically the full infinite and continuous set of GSs configurations by means of displacements at no energy cost. This zero-frequency freely-sliding mode is the Goldstone mode consistent with an emerging continuous translational invariance of the model, connecting continuously with an acoustical phase mode (phason) at finite wavelength. By contrast, in the pinned soft-chain region K < K c , all particles remain trapped close to the substrate-potential minima, a configuration which exhibits a finite energy barrier against motion over the corrugation.
The locking is provided here, despite translational invariance, by the inaccessibility of forbidden configurations, which act as dynamical constraints. Above K c the incommensurate-chain sliding can therefore be initiated by an arbitrarily small driving force, whereas for K < K c the chain and the corrugation lock together through the pinning of the kinks or superkinks separating locally lattice-matched regions. Note that the locking of the free ends in a finite-size FK chain necessarily leads to pinning, even when θ is irrational and regardless of how large K is. Even for a finite chain, it is nevertheless still possible to define and detect a symmetry-breaking Aubry-like transition [32][33][34].
For the characterization of the Aubry transition, a "disorder" parameter ψ was conveniently defined [35] as the smallest distance of atoms from the nearest maximum of the corrugation potential. This quantity vanishes in the freely-sliding state, and is nonzero in the pinned state. At the critical pinned-to-sliding point, the disorder parameter exhibits a power-law behavior Here the values of the critical exponents are functions of the irrational length ratio θ. Specifically, for the golden-mean ratio χ ψ ≃ 0.7120835 and χ PN ≃ 3.0117222 [31,[35][36][37][38][39][40][41]. Equation (15) characterizes the continuous Aubry transition with a scaling behavior typical of critical phenomena, here at T = 0 but as a function of the the stiffness parameter K. It is common to refer to the exponents in Eq. (15) as super-critical, since they are specific to the pinned side of the transition, K ≤ K c . Sub-critical exponents were introduced for the freely-sliding state K > K c , as well. To describe the response of the model to an infinitesimally small dc force F applied to all atoms, an extra damping term −mγẋ i has to be included in the equation of motion to prevent unlimited acceleration, and to achieve instead a steady-state. The resulting effective viscosity in the subcritical region is defined as Γ = lim F →0 F/(mv), in terms of steady-state average velocity v resulting in response to F . At the Aubry critical point K c , the effective viscosity Γ diverges. For θ = (1 + 5 1/2 )/2, the golden-mean ratio, the scaling behavior of Γ is with χ Γ ≃ 0.029500. As is the case for all scaling relations, Eq. (16) provides the leading divergence close to the Aubry point; at a larger distance from K c , Γ deviates from Eq. (16). Eventually, in the K → ∞ SG limit, Γ decreases toward γ. In general, in the unpinned phase at K > K c the incommensurate FK model exhibits an effective viscosity systematically larger and thus a mobility B = v/F = (mΓ) −1 consistently smaller than its maximum value (mγ) −1 . This observation is illustrated by the F → 0 limiting values of the K > K c curves of Fig. 3. Exclusively in the K → ∞ SG limit, the incommensurate system moves under an infinitesimal force without any extra dissipation added to the base value γ, therefore in a frictionless sliding motion, despite the finite corrugation magnitude U 0 . The first prediction of vanishing static friction was formulated for the incommensurate infinite-size sufficiently hard FK chain by Peyrard and Aubry [31]. This phenomenon was subsequently re-discovered for incommensurate tribo-contacts, and named superlubricity [42,43]. This name has drawn criticism, because it could misleadingly suggest the vanishing of kinetic friction too, in analogy to superfluidity or superconductivity. Actually instead, the depinning of sliding interfaces closes just one of the channels for energy dissipation, namely the one associated to the stick-slip elastic instability at low speed. Additional dissipation channels, including the emission of vibrations such as sound waves, remain active, with the result that the actual kinetic friction force remains nonzero and growing with increasing sliding speed. All the same, the superlubric state does attain a significant reduction of the kinetic friction force (and thus an increased mobility B), compared to the pinned state K < K c , see Fig. 3. The driven FK model was usefully employed to describe the onset of sliding of a crystalline contact [44], even though this model cannot describe rigorously the experimentally-significant real-life plastic deformations of the contact.
Experimentally, superlubricity has been studied for a graphite flake sticking to the tip of an atomic force microscope (AFM) sliding over an atomically flat graphite surface [45][46][47]. Extremely weak friction forces of less than 50 pN were detected in the vast majority of the relative flake-substrate orientations, namely those orientations generating incommensurate contacting surfaces, see Fig. 4. Stick-slip motion, associated with a much higher friction force (typically 250 pN), was instead found in the narrow ranges of orientation angles where the flake-substrate contact was commensurate.
The above discussion ignores temperature, assuming so far T = 0. At nonzero temperature, the sliding-friction response of the FK model [48,49] requires of course the addition of a thermostat (see Sec. 2). The common choice of a Langevin thermostat for example simulates all dissipation mechanisms through a viscous force −mγẋ i , and includes fluctuations by the addition of Gaussian random forces whose variance is proportional to temperature T , as sketched in Sec. 2.1. At T > 0, thermal fluctuations can always overcome all sorts of pinning and will thus initiate sliding by nucleation of mobile defects even in the fully commensurate (pinned) condition, see Fig. 2.
More generally, in the FK model the dimensionless coverage θ = a b /a c plays a central role, because it defines the concentration of "geometrical kinks" close to θ = 1 and of "superkinks" which arise when θ deviates slightly from a background commensurate pattern that is not θ = 1, but rather a rational θ 0 = p/q, with p and q mutually prime integers. If θ is only slightly different from θ 0 , the GS of the FK model consists of extended domains with the commensurate pattern associated to θ 0 , separated by superkinks (super-antikinks), in the form locally mismatched regions of compression (expansion) relative to θ 0 .
The above concepts of pinning or superlubricity apply for an infinitesimal applied force. Additional interesting physics arises at finite force. By increasing the driving force, a FK model with a pinned GS (either commensurate, or incommensurate but past the Aubry transition) is known to show a hierarchy of first-order dynamical phase transitions, starting from the completely immobile state, passing through several intermediate stages characterized by different running states of the kinks, to eventually reach a totally running state. Consider, for example, the ratio θ = 21/41: initially the mass transport along the chain is supported by superkinks constructed on top of the background θ 0 = 1/2 structure. Since the average superkink-superkink distance 41 a b is large, they interact weakly, and the atomic flow is restricted by the need for these rarefied superkinks to negotiate their PN barriers (see Fig. 5a). For larger driving F , the effective PN barriers are tilted and lowered (in analogy to barriers of the corrugation potential), producing an increased single-kink mobility B. As a result, the zero-temperature transition from the locked state (B = 0) to the running state takes place at the force F ≈ F tk = Cπε PN /a b , where the factor C ∼ 1 depends on the shape of the PN potential. In terms of the dimensionless superkink concentration θ k = 1/41, the mobility becomes B ≃ θ k B f .
Beyond F tk , further possibilities depend on the damping coefficient γ. At very small damping, γ < 0.05, the driven model transition leads directly into the fully running state, because running superkinks self-destroy soon after they start to move, causing an avalanche, thus driving the whole chain to a total running state similar to that shown at the right side of Fig. 2. When the dissipation rate is larger, γ > 0.05, one instead finds intermediate stages with stable running superkinks, see Fig. 5. The mechanism for a second rapid increase of the mobility after depinning depends again on the value of γ (for details see Refs. [15,50,51]). In between the initial superkink-sliding stage and the fully running state, a sort of "traffic-jam" intermediate regime may emerge [15].
A qualitatively similar picture was confirmed also for different and more complex kink patterns, such as that shown for the θ = 21/31 example in Fig. 5b [50]. In this case, the GS consists of domains of the θ 0 = 2/3 commensurate structure, separated by superkinks at an average spacing 30 a b . Even the θ = 2/3 pattern could itself be viewed as a dense array of trivial kinks constructed on top of the simple θ 0 = 1/2 background structure. The force dependence of the mobility B(F ) bears a trace of this double nature, with a state of running superkinks preceding a state of running kinks.
For not too small γ therefore the mobility B increases in two distinct steps as the driving force is increased. A first step, at F = F sk ≈ 0.08 F 0 (here F 0 = πU 0 /a b defines the depinning force for the fully commensurate model θ = 1) occurs when the superkinks begin to slide; then a second step, at F = F ′ tk ≈ 0.18 F 0 , occurs in correspondence to the unpinning of the trivial kinks.

Extensions of the FK Model
Several extensions of the FK model have been proposed to describe a broad range of frictionally relevant phenomena. Most of these generalizations involve modifications of either the interactions or the system dimensionality. To address more realistic systems, anharmonicity of the chain interatomic potential has been studied in detail. The resulting features include mainly new types of dynamical solitons (supersonic waves), a modification of the kink-kink interaction, the breaking of the kink-antikink symmetry, and even the possibility of a chain rupture associated to the excessive stretching of an antikink [15]. The large kink-antikink asymmetry consistent with friction experiments in layers of repulsive colloids [52] was attributed to the strong anharmonicity in the colloid-colloid interaction [53,54]. The essence of this asymmetry is the same as that between the physical parameters of a vacancy and those of an interstitial.
Research has addressed also substrates with a complex corrugation pattern [22,55], including quasiperiodic [56,57] and random/disordered corrugation profiles [58][59][60]. Modifications from the plain FK model may generate qualitatively different types of excitations, e.g. phonon branches and kinks of different kinds, as well as modifications in the kink-antikink collisions. At small driving force, where the dynamics and tribology are dominated by moving kink-like structures, different sliding modes appear.
The Frenkel-Kontorova-Tomlinson (FKT) model [61,62] introduces an harmonic coupling of the sliding atomic chain to a driving support, thus making it possible to investigate stick-slip features in a 1D extended simplified contact. The FKT framework provided the ideal platform to investigate the tribological consequences of combined interface incommensurability, finite-size effects, mechanical stiffness of the contacting materials, and normal-load variations [63,64].
Important generalizations involving increased dimensionality compared to the regular FK model bear significant implications for tribological properties such as critical exponents, size-scaling of the friction force, depinning mechanisms, and others. In particular 2D extensions of the FK model [4,15] have been applied to the modeling of the (unlubricated) contact of two crystals. Such is the case, for example, in quartz-crystal microbalance (QCM) experiments, where single-layer adsorbate islands are made slide over a crystalline substrate [65]. Another case is that of recent experiments carried out with 2D monolayers of colloids driven over a laser-generated optical lattice [52][53][54][66][67][68][69][70].
Interestingly, the 2D Aubry transition of incommensurate colloids was shown by Mandelli et al. [70] to be of first order, rather than of second order as in 1D. As a consequence in 2D the free sliding and the pinned phases retain local stability for a range of parameters that extends beyond the transition point, a point where the total energy has a crossing singularity instead of a smooth stiffness dependence as in 1D. It is likely, although not proven to our knowledge, that the 2D FK should possess a first-order Aubry transition too.
These approaches which generalize the FK model have been of use for the study of the transient dynamics at the onset of sliding. Capturing these transient phenomena is often highly nontrivial in fully realistic MD simulations (see e.g. Ref. [83]). An interesting example of such a transient is the depinning of an atomic monolayer driven across a 2D periodic substrate profile of hexagonal symmetry [83]. The formation, by nucleation, of an island of moving atoms in a "sea" of quasi-stationary particles mediates the transition from the locked to the running state. The moving island expands rapidly along the direction of the driving force, and grows at a slower rate in the orthogonal direction. Within the island, the 2D crystal retains its approximate ordered hexagonal structure, thanks to its stiffness supported by the intra-atomic forces. As a result, at the onset of depinning, the model exhibits regions of almost perfect hexagonal-lattice order delimited by a closed boundary of dislocations.
Kinks in 1D and dislocation lines in 2D exhibit peculiar tribological properties. In noncontact experiments, an oscillating AFM tip was seen to dissipate significantly more when hovering above a dislocation line of a incommensurate adsorbates than above inregistry regions [84]. The larger softness and mobility of the dislocation regions [85][86][87][88] accounts for this effect. An explicit demonstration of this mechanism was carried out by the study of a mismatched FK chain, whose dynamics was forced and simultaneously probed by an oscillating localized model tip [89]. This approach illustrates the ability of the FK model to capture the local modifications of the dissipation properties. In contrast, if retardation effects related to the finite speed of sound across a material need to be taken into account, more sophisticated models are called for.

Quantized sliding velocity
The investigation of systems confined between two shearing sliders, such as single particles [90][91][92] or harmonic chains [93][94][95] embedded between competing periodic potentials, have led to the discovery of several nonlinear tribological phenomena involving either stick-slip dynamics or the formation of peculiar "synchronized" sliding regimes [96,97]. The FK model can be generalized with the addition of a second, different sinusoidal corrugation potential, as sketched in Fig. 6c. When the second potential is spatially advanced relatively to the first as a function of time, the model realizes the simplest idealization of a slider-solid lubricant-slider confined geometry. In this extended model the lattice mismatch was shown to generate peculiar and robust "quantized" sliding regimes [95,[98][99][100][101][102][103][104][105], where the chain deformations are synchronized to the relative motion of the two corrugations, in such a way that the chain's (i.e. the solid lubricant's) average drift velocity acquires nontrivial fixed ratios to the externallyimposed sliding velocity. Specifically, the ratio of the lubricant speed to that of the slider w = v cm /v ext remains locked to specific "plateau" values across broad ranges of most model parameters, including the potential magnitude of the two sliders, the chain stiffness (see Fig. 6a), the dissipation rate γ, and even the external velocity v ext itself. The speed ratio is ultimately determined by geometry alone: w = 1 − r −1 + , where r + = a + /a 0 is the incommensurability ratio between the chain spacing a 0 and that, a + , of the closest slider [98]. The plateau mechanism is the fact that the kinks formed by the mismatch of the chain with one slider (the slider whose spatial periodicity is closest to that of the chain), are rigidly dragged at velocity v ext by the other slider. The kink density being geometrically determined and lower than the chain density implies that the overall velocity ratio shares exactly the same properties. The exactness of the velocity plateaus implies a sort of "dynamical incompressibility", identically zero compliance to perturbations trying to modify v cm from its quantized-plateau value. This robustness of the plateaus can be demonstrated e.g. by adding a constant force F ext , pushing all particles in the chain: as long as F ext is small enough, it does perturb the dynamics of the velocity-plateau attractor, but not the value of v cm . Eventually, above a critical force F c , the driven model abandons the plateau dynamics. This transition, explored by increasing the external driving force F ext , exhibits a broad hysteresis, and shares many features of the static-friction depinning transition, except that here it takes place "on the fly" [101,103]. Disregarding details, this transition is then formally equivalent to the standard Aubry depinning transition [30,31], with the moving kinks of the lubricantsubstrate interface taking here the role of particles. The robustness of the quantized plateau stands even after replacement of the sinusoidal corrugation potential of Eq. (14) with a deformed profile: the Remoissenet-Peyrard non-sinusoidal potential even extends the velocity plateau in the space of model parameters [106].
The quantized sliding regime of the crystalline solid lubricant was also investigated in a significantly less idealized 2D model, including the perpendicular degree of freedom [107]. MD simulations carried out for a monolayer or multilayer lubricant film where atoms interacting via Lennard-Jones potentials can also move perpendicularly to the sliding direction, as sketched in Fig. 7 showed quantized plateaus in this case too. These plateaus were shown to be resilient against a variations in the loading forces across a broad range, against thermal fluctuations, and also against the presence of quenched disorder in the substrates. This quantized sliding state was also characterized by significantly lower values of kinetic friction F k (the average force needed to maintain the advancement of the top slider) [108], than the regular non-quantized regime, see Fig. 8.
Quantized sliding has been again demonstrated more recently in a 3D model where the lubricant is represented by a layer of Lennard-Jones atoms [109]. The quantized-sliding state and its boundaries were fully characterized in the special case of perfectly-aligned crystalline layers. We note however there are reasons to expect that incommensurately mismatched epitaxial layers could relax to a mutually rotated alignment [110,111]. The quantized-sliding state in such rotated arrangements is the subject of active investigation. More generally, no experimental observations of the quantized sliding predicted for solid lubricants has appeared so far. Layered systems such as graphene and BN appear to offer a good opportunity for the future study of these curious phenomena.

Molecular Dynamics Simulations
The simple models considered above yielded precious qualitative and often semiquantitative understanding of several features of friction. To address subtler physical behavior in specific systems, it is nevertheless necessary and desirable to include atomistic structural and mechanical details of the interface. MD simulations can help make an inroad such detail, also offering a level of detail that can in some instances replace experiment.
Thanks to advances in computing algorithms and hardware, recent years have witnessed a remarkable increase in our ability to simulate tribologic processes in realistic nano-frictional systems, and obtain detailed microscopic information. A MD simulation is de facto a controlled computational experiment, where the overall atomic Figure 9. Two sketches of typical MD simulations. Left: the melting induced by shearing at a boundary-lubricated interface [112]. Right: simulations of a monolayer of colloids interacting with each other and with a periodic corrugated potential produced by a laser interference pattern as in Ref. [53].
dynamics is provided by the numerical solutions of suitably generalized Newton's equation of motion relying on interatomic forces derived by specific realistic interparticleinteraction potentials. Tribological simulations require a careful selection of the geometric arrangement of the sliding interface, e.g. as in Fig. 9, and of the applied boundary conditions. Influential review articles [10,11] cover the atomistic MD simulation of friction, with a focus on technical aspects such as the construction of a realistic interface, suitable techniques for the application of load, shear, and temperature control. The simplest approach to temperature control, namely adding a Langevin thermostat to Newton's equations [9], has been adopted broadly, but more refined approaches have been proposed and adopted for friction simulations, as discussed in Sect. 4.1 below. Physically relevant quantities, including the average friction force, the slider and the lubricant mean velocities, several correlation functions, and the heat flow can be evaluated numerically by carrying out suitable averages over the model dynamics of a sliding interface, as long as it is followed for a sufficiently long time. The modeling of friction must first of all address correctly ordinary equilibrium and near-equilibrium phenomena, where the fluctuation-dissipation theorem (Sec. 2) governs the smooth conversion of mechanical energy into heat, but most importantly it must also deal with inherently nonlinear dissipative phenomena such as instabilities, stick-slip, and all kinds of hysteretic response to external driving forces, characteristic of non-equilibrium dynamics.
The choice of realistic interatomic forces is often a major problem. Indicating with U{R 1 , R 2 , ...R N } the total interaction energy as a function of all atomic coordinates {R i }, the force on atom i is F i = −∇ R i U, fully determined in terms of U. Unfortunately, the adiabatic energy U results from the solution of the quantum ground state of the electrons -a practically complicated problem whose quantitative outcome may moreover be of uncertain quality. The reason why ab-initio MD, e.g. of the Car-Parrinello type [113], is not generally used in sliding friction is that it can neither handle large systems, exceeding few hundreds of atoms, nor run for tribologically significant duration, usually in excess of ∼ 1 ns. On the other hand, the physical situations where a first-principles description of interatomic forces is mandatory are not too common in the frictional phenomena studied so far. As a consequence, most MD models for friction rely on more or less refined interatomic "force fields", ranging from sophisticated energy surfaces modeled on calculations at ab-initio density-functional or tight-binding level [114,115], to empirical distance-and angle-dependent manybody classical potentials, to basic pairwise potentials (e.g. Morse or Lennard-Jones), to the simplest elastic-springs models, which represent generalizations of the FK model. Concretely, the scientific literature documents many realistic force fields, ready to address several classes of materials and their combinations [116,117]. While these force fields allow qualitative atomistic simulations of tribological systems, their limitations often prevent quantitative accuracy. In particular, in the course of such a violent frictional process as wear, atoms are likely to modify their chemical coordination and even their charge state: phenomena and radical chemical changes usually impossible to describe with empirical force fields. Mechanochemistry and tribochemistry are timehonored areas [118,119] offering obvious examples where empirical force-fields would fail, and simulations must by necessity be conducted by electronic-structure based first-principles methods. Also, even if for a specific system, element or compound, a satisfactory force field has been arrived at, the mere replacement of one atomic species with another one generally requires a complete and usually difficult re-parameterization of the whole force field. As a result, quantitatively accurate nanofrictional investigations remain a substantial challenge because of opposite limitations in the use of first principles versus empirical force fields. A promising compromise could possibly be provided by the so-called reactive potentials [120][121][122], capable of describing some chemical reactions, including interface wear with satisfactory computational efficiency in large-scale atomic simulations, compared to semi-empirical and first-principles approaches.
Retardation effects due to the finiteness of the speed of sound are usually irrelevant in slow-speed experiments (v < 1 mm/s). For larger speed, retardation effects related to the finite speed of sound can be taken into account explicitly in MD modeling, provided that rigid layers are either omitted or introduced with special care. Other effects of nonlocality in time, such as retardation due to the finite speed of light [123,124] are usually omitted in the MD force fields altogether, as they lead to negligible corrections in all conditions where sliding involves a proper material contact. Such retardation effects do play a role in noncontact geometries, such as in experiments probing lateral Casimir forces [125], whose strength can become relevant at large sliding speed [126].

Thermostats: the Dissipation of Joule Heat
As we already mentioned above in Sec. 2, any kind of sliding friction involves mechanical work, some of which is then transformed into heat (the rest going into structural transformations, wear, etc.). The heat is then transported away by phonons (and electrons in the case of metallic sliders) and eventually dissipated to the environment. Likewise, all excitations generated at the sliding interface in simulations should be allowed to propagate away from it, and to disperse in the bulk of both sliders. Instead, due to the small simulation size, this energy may unphysically pile up in the rather small portion of solid representative of the "bulk" of the substrates, where these excitations are scattered and back-reflected by the simulation-cell boundary, instead of being properly dissipated away. In order to prevent continuous heating and attain a steady state of the tribological system, the Joule heat must then be removed at a steady rate. In the FK and PT models, a viscous damping term −mγẋ i , Eq. (3), is generally introduced for this purpose. In these minimal models however, the value of γ is well known to affect the dynamical and frictional properties, but there is unfortunately no clear prescription for the choice of γ. In MD atomistic simulations, the heat removal is often achieved by means of equilibrium thermostats, e.g. Nosé-Hoover or Langevin, see Sec. 2.1. In this way however an unphysical energy sink is spread throughout the simulation cell. As a result, the atoms at the interface fail to follow their actual conservative trajectories, but evolve through an unphysically damped dynamics, with unknown and generally undesired effects on the overall tribological properties [127]. In order to address and mitigate this problem, modifications of the equations of motion for the atoms inside the microscopically small simulation cell were proposed with the target of reproducing the frictional dynamics of a realistic macroscopic system, after the integration of extra "environment" variables. One possible approach is the application of Langevin equations with a damping coefficient that changes as a function of the position and velocity of each atom in the lubricant, in accordance with the dissipation known for the atoms adsorbed on a surface [128]. This method involves modifying the standard Langevin equations [129]. Another approach to improve the simulation of dissipation within blocks in reciprocal motion requires modifying the damping term (3) to a form where v loc is the average center-mass velocity of the atoms forming the sliding block to which particle j belongs locally [130,131]. Another more rigorous, physically appealing approach is the recently-implemented dissipation scheme, drawing on earlier, longknown formulations [132][133][134] and subsequent derivations [135][136][137], describing the correct embedding of the Newtonian simulation cluster inside a larger heat bath made of the same material. Upon integrating out the heat bath degrees of freedom, atoms in the boundary layer that borders between the cluster and the heat bath are subjected to additional non-conservative and non-Markovian forces that mimic the surrounding bath through a so-called memory kernel. An approximate but very practical scheme replaces this memory kernel by a simple viscous damping γ, here applied exclusively to atoms in the boundary layer. The magnitude of the parameter γ is optimized variationally by minimizing, with surprising accuracy, the energy reflected across the boundary [138,139]. This dissipation scheme has been implemented recently in nanofriction simulations where it was shown to improve greatly over other conceptually and practically inadequate thermostats. Besides the limitations of system size and simulation times that are obvious and will be discussed later, there is another limitation concerning temperature, that is rarely mentioned. All classical frictional simulations, atomistic or otherwise, are only valid at sufficiently high temperature. They become in principle invalid at low temperatures where the mechanical degrees of freedom of solids progressively undergo "quantum freezing", and both mechanics and thermodynamics deviate from classical. Unfortunately there is at the moment no available route to include appropriately these quantum effects in dynamical and frictional simulations.

Size and timescale issues
Each core of a present-day CPU executes ∼ 10 9 floating-point operations per second (FLOPS). MD simulations usually benefit effectively of medium-scale parallelization. Approximately linear scaling can be achieved up to ∼ 100 cores, thus a MD simulation can execute ∼ 10 11 FLOPS routinely. The evaluation of the forces is usually the most CPU-intensive part of a MD simulation. For each atom, depending on the force-field complexity and range, this evaluation can require ∼ 10 − 10 2 operations, or even more. As a result, the number of time-integration steps N step multiplied by the number of simulated particles N, is ∼ NN step ≃ 10 10 per computer runtime second. Given that simulations of atomic-scale friction require time-steps in the femtosecond region, a medium-size simulation involving N = 10 5 simulated particles, can advance at an estimated speed of ∼ 10 5 fs each real-life second, namely ∼ 10 9 fs = 1 µs each simulation day. Clearly, speed scales down for more refined force fields, and for larger systems size, although this increase may be mitigated by a larger-scale parallelization.
We can compare these estimates with typical sizes, duration, and speeds in frictional experiments. In macroscopic tribology experiments, sliding speeds often range in the 0.1 − 10 m/s region: each microsecond the slider would progress by 0.1 to 10 µm, namely ∼ 10 3 − 10 4 lattice spacings of standard crystalline surfaces. In such conditions 10 −3 − 10 −2 µs may suffice to generate a good statistics of atomic-scale events, although it may still be insufficient to address e.g. the diffusion of wear particles or additives in the interface, or phenomena associated to surface steps and/or point defects. By contrast, in nanoscale AFM experiments the tip usually advances at much lower speeds ≃ 1 µm/s: over a typical run it is possible to simulate a tiny ∼ 1 pm displacement, far too small to explore even a single atomic-scale event, let alone averaging over a steady state. For this reason, in all conditions where long equilibration times and/or slow diffusive phenomena and/or long-distance correlations can be expected, models should be preferred to realistic but expensive MD. However, MD simulations can provide so much physical insight that they make sense even if carried out at much higher speeds than in real-life AFM or surface force apparatus (SFA) experiments: in practice, currently the sliding speeds of most atomistic tribology simulations are in the ∼ 1 m/s region.
Here however we should distinguish between static and kinetic friction, and for the latter between smooth-sliding and stick-slip regimes. Smooth kinetic friction generally increases with speed (velocity strengthening), but sometimes decreases with increasing speed in certain intervals. In the former case, simulating smooth high-speed frictional sliding is not fundamentally different from the real sliding at low speed, with appropriate changes in frictional forces with v. Velocity weakening conditions, alternatively, tend to lead to an intrinsic instability of smooth sliding, which is therefore not often pertinent to real situations. As a result, for nanoscale systems, MD simulations is of value in the description of smooth dry kinetic friction despite the huge velocity gap. On the other hand, static friction -the smallest force needed to set a slider in motion -is also dependent on the simulation time (a longer wait may lead to depinning when a short wait might not), and generally dependent on system size, often increasing with sub-linear scaling with the slider's contact area. To address this kind of behavior in MD simulations, it is often necessary to resort to scaling arguments in order to extrapolate the large-area static friction from small-size MD simulations [131,140].
Returning to the simulation time problem, let us come to stick-slip in MD simulation, and to the desirability to describe the stick-slip to smooth-sliding transition as a function of parameters such as speed. In AFM and SFA experiments, stickslip and its associated characteristically high friction and mechanical hysteresis tend to transition into smooth sliding when the speed exceeds ∼ 1 µm/s; in contrast, in MD modeling the same transition is observed in the ∼ 1 m/s region. This 6-order-ofmagnitude discrepancy in speed between experiments and simulations is well known and has been largely discussed [141][142][143][144] in connection with the effective mass distributions and spring-force constants, that are vastly different, and highly simplified in simulations. Attempts to fill the time and speed gaps can rely on methods, such as hyperdynamics, parallel-replica dynamics, on-the-fly kinetic Monte Carlo, and temperature-accelerated dynamics which have been developed in the last decades [145][146][147]. However, caution should generally be exerted in that some of these schemes and methods are meant to accelerate the establishment of equilibrium but not always to treat properly the actual frictional-loss mechanisms. Concerning stick-slip friction, another problem is that, unlike simulations, real experiments contain mesoscale or macroscale component intrinsically involved in the mechanical instabilities of which stick-slip consists. Here the comforting observation is that stick-slip is nearly independent of speed, so that so long as a simulation is long enough to realize a sufficient number of slip events, the results may already be good enough [148]. One can even describe stick-slip friction adiabatically, e.g., from a sequence of totally static calculations, where a periodic backand-forth sliding path is trodden, the area of hysteresis cycle generated by two different to and from instabilities representing the friction [149].
A serious aspect of stick-slip friction which MD simulation is unable to attack is ageing. The slip is a fast event, well described by MD, but sticking is a long waiting time, during which the frictional contact settles very slowly. The longer the sticking time, the larger the static friction force necessary to cause the slip. Typically experiments show a logarithmic increase of static friction with time [150]. Rate and state friction approaches, widely used in geophysics [151], describe phenomenologically frictional ageing, but a quantitative microscopic description is still lacking. Mechanisms invoked to account for contact ageing include chemical strengthening at the interface in nanoscale systems [152], and plastic creep phenomena in macroscopic systems [153]. Contact ageing is observed also in other disordered systems out of equilibrium, including glasses and granular matter. In seismology finally, as will be discussed later, it is generally accepted that ageing is responsible for aftershocks, as also shown by some models [154].

Multiscale Models
If MD simulation may be satisfactory in nanoscale friction, it is clearly not capable of describing mesoscale and macroscale tribology. The insurmountable difficulties of the fully atomistic treatment of all typical and large length scales that are responsible for the dynamical processes in large scale systems has prompted in recent years increasing efforts towards multiscale approaches. The main idea is that, at a sufficient distance from the sliding interface, continuum mechanics should describe all processes to a fair level of approximation. Finite-elements simulations of the continuum mechanics may provide a practical model for the elastic and plastic deformations. Using finite element methods, one can increase the coarse-graining level while moving away from the sliding interface, thus keeping the computational effort under control. Several groups [155,156] combined the MD description of the sliding interface, where local deformations at the atomic length scale and highly nonlinear phenomena occur, with a continuum-mechanics description in the "bulk" regions where strains are continuous and small. The main difficulty faced by this class of approaches is the correct choice of the matching between the atomistic region with the continuum part [157]. Because at the continuum level the detail of lattice vibrations cannot be represented in full, the matching conditions should at least minimize the reflection of the acoustic phonons at the atomistic-continuum interface. In other words, the matching should allow the transmission of sound deformations in both directions with sufficient accuracy: this is necessary for a proper disposal of the Joule heat into the bulk.

Selected Results of MD Simulations
Simulations can provide direct insight in the dynamical processes at the atomistic level, that are at the origin of friction, allowing a connection of these microscopic facts with their macroscopic counterparts. Case studies, in which the system is well describable by both experimental and theoretical sides, are of extreme importance firstly to permit a crosscheck between the two, and then to make use of simulations in order to highlight particular aspects that cannot be accessed by experiments. Here we summarize the main results of a few selected simulations sampled from the expanding literature of friction simulation, certainly not claiming an exhaustive review of the field.

Rare-gas islands and metal clusters.
The sliding of rare-gas overlayers deposited on metallic substrates at low temperature has contributed much to the understanding of how friction scales with the contact-area size, the substrate corrugation, and the sliding velocity. Rare-gas atoms condense into 2D solid islands showing a facetedcircular shape, arranged on multiple layers at low temperatures or on a single layer at diffusion-enabled temperatures [158]. The friction characteristics of these solid islands on the substrate, resulting from inertial sliding, has been probed experimentally by QCM apparatuses, revealing a complex interplay among friction, coverage, and temperature [131,159]. The rare-gas lattice spacing inside the island, generally incommensurate but sometimes commensurate with that of the substrate, plays a very important role in determining the pinning or free sliding that controls the frictional behavior. A generally overlooked aspect which has been highlighted only recently [131] is the larger thermal expansion coefficient of rare-gas layers than that of a metal substrate, causing a temperature-dependent lattice mismatch at the interface with possible incommensuratecommensurate transitions. Due to this mechanism, MD simulations have predicted the possible appearance of static-friction peaks in the correspondence of a long-range commensurate phase occurring at a particular temperature [24]. In the case of periodic monolayers, a change in the lattice mismatch can be also induced by an adhesion-driven densification of the adsorbate [160], again eventually encountering a commensurate phase [131,159]. Simulations of rare-gas incommensurate adsorbates, whose linear substrate-induced misfit dislocations ("solitons") must flow during sliding, have revealed the role of their entrance in the depinning of the island [24], and of their dissipation through anharmonic coupling to phonons, in kinetic friction [160]. Finite-size effect are in this case of absolute relevance, effectively generating (or enhancing) static friction through a pinning barrier arising at the interface edge, which solitons must overcome to establish motion [24]. The edge-related origin of the pinning mechanisms implies that static friction F s can grow with the island size A at most as F s ∝ A 1/2 , i.e. as its perimeter P ∝ A 1/2 , if the pinning points were uniformly distributed along the island or cluster edge. As shown in Fig. 10, a different shape of the deposited nano-object can generally lead to a different scaling exponent. Similar sublinear scaling exponents were identified in dynamic friction experiments in which gold nanoclusters of variable size/shape were dragged at low speed over a graphite substrate by an AFM tip [161]. Scaling exponents of both the rare-gas island/metal surfaces (theoretical) and dragged gold clusters/graphite (experimental) are in the order of 1/4. This indicates that not all points at the boundary provide pinning with equal efficiency. A scaling close to ∼ A 1/4 might rather indicate a random efficiency of boundary points, whereby only ∼ P 1/2 provide effective pinning.

AFM, nanotubes, and other systems.
Nowadays' computational capabilities even permit the atomistic simulation of an entire AFM tip, enabling the understanding of several mechanisms which are not describable by simplified PT-like models (see Fig. 11). For example, it is possible to highlight the formation/rupture dynamics of contacts in multi-asperity interfaces, and consequently estimate the true contact area as a function of the apparent one. Besides, it is possible to investigate the effect of the tip plasticity and elasticity, which are of fundamental importance to define the load-dependent contact area [162], and as channels for dissipation and wear [163][164][165][166][167]. This approach enables the bottom-up derivation of the linear scaling laws of macroscopic friction with size, and their transition to the sublinear ones for incommensurate nanosized contacts. We can now understand that such transition takes place when the contact roughness becomes large compared to the range of interfacial interactions [162].
In the study of repeated scratching of metallic surfaces by hard AFM tips, widely employed in the field of micro/nano machining, MD simulations have uncovered strongly non-linear trends of the frictional force with the feed (i.e. the distance from the first groove), induced by lateral forces exerting on the tip due to the substrate plasticity [168].
It is also important to mention the simulations of nanotubes (NT), either made of carbon or hexagonal BN, which, due to their extraordinary mechanical and electronic properties, have been investigated with enormous interest in the last decades. Almost defectless NT can be formed nowadays with lengths of the order of 1 cm [169], and precise measurements of their mechanical and frictional properties have started to appear in literature [170,171]. Simulations of concentric nanotubes in relative motion (telescopic sliding), have revealed the occurrence of well-defined velocities at which friction is enhanced, corresponding to a washboard frequency resonating with longitudinal [172] or circular [173] phonon modes, leading to enhanced energy dissipation. The frictional response becomes highly non-linear while approaching the critical velocity and, contrary to macroscopic systems, washboard resonances can arise at multiple velocities, especially for incommensurate interfaces where more than one length scale may be in common to the contacting surfaces [172].
The exceptional electro-mechanical properties of NTs have also been investigated by various tip-based techniques, revealing a strong friction anisotropy dictated by NTs orientation. In this respect, simulation-assisted experiments of a sliding nanosized tip over CNTs reveal that transversal friction is enhanced by a hindered-rolling motion of the NT, with a consequent frictional dissipation that is absent in the longitudinal sliding [174]. The same elastic deformation have been reportedly responsible of a reverse stick-slip effect in the case of an AFM probe sliding over a super-lattice CNT forest [175]. Here, simulations reveal that the fast sticking is induced by the penetration of the tip into the valley between the NTs and its interaction with the tubes on both sides, causing an elastic shell buckling of the CNTs. In contrast, the gradual slipping occurs over much longer distance because it includes both the sliding on the top of the NT and the energy release at both sides of the graphitic wall.

4.4.3.
Boundary lubricated sliding. When two sliding surfaces are separated by a thick lubricant film, as it ordinarily happens under weak-load conditions, the tribological response of the confined system is typically determined by the fluid viscosity. In these cases of hydrodynamic lubrication, friction can be computed based on the Navier-Stokes equations, which prescribe a monotonically increasing kinetic friction as a function of the relative sliding speed [176]. By contrast, at high load and low driving velocity, the lubricant may not maintain a broad gap between the sliding surfaces, with the result that solid-solid contact eventually occurs. Prior to full squeezeout under pressure, as confirmed by experiments and simulations, the intervening boundary film usually changes from liquid to solid or nearly solid, exhibiting a layered structure prone to develop finite static friction and a high-dissipative stick-slip dynamics in such "boundary-lubrication" regime. Both SFA measurements [177] and MD investigations [178,179] have demonstrated sharp upward jumps of friction at squeezeout transitions where the number of lubricant layers decreases from N to N − 1. Boundary-lubricated systems often display stick-slip dynamics during tribological measurements associated with a significant value of the friction dissipation. As the load increases, however, it becomes harder and harder to squeeze out an extra lubricant layer. This hardening and increased difficulty of squeezeout reflects the tendency to crystallization of the initially liquid lubricant and the increased cost of the "crater" whose formation constitutes the nucleation barrier of the transition [180][181][182]. Once it has happened, generally this N → N − 1 re-layering transition gives rise to an upward friction jump. In principle however, upon re-layering of the solidified (structured) confined film, the twodimensional (2D) parallel crystalline-like order could occasionally change under pressure toward a more favorable mismatched (incommensurate) substrate-lubricant geometry. In that case, as a result, sliding friction might actually switch downward from highlyfrictional stick-slip to smooth dynamical regimes characteristic of incommensurate superlubric interfaces, with a highly mobile 2D soliton pattern, of the type sketched in Fig. 12. So far this type of event has only been observed in simulations [112].
In SFA experiments, boundary-lubricated systems often display stick-slip dynamics during tribological measurements, associated with a significant value of the friction dissipation. One cannot directly access the detailed film and interface rearrangements giving rise to the stick-slip mesoscopic intermittent dynamics. The mechanisms at play for the stick-slip dynamics in the boundary-lubrication regime have been studied by MD investigations [183][184][185][186]. Several realistic models for lubrication layers were simulated [144,[187][188][189][190]. The issue whether frictional shearing occurs through the middle of the solid lubricant film, possibly accompanied by melting-freezing, or if it forms a smooth shear band, or it occurs at the substrate-lubricant boundary, is one which can be addressed by computer modeling. Depending on the relative strength of the potentials governing the lubricant-lubricant and lubricant-substrate interactions, a thin confined film may exhibit a solid-like or liquid-like behavior under shear. If the interaction with the substrates is weaker than the lubricant-lubricant one, then sliding takes place mainly at the surface-lubricant interfaces. The lubricant film is then allowed to maintain or acquire a solid order. If both the solid lubricant and the substrate are characterized by nearly-perfect crystalline structures and these structures are mismatched and/or misaligned [111,[191][192][193][194], then smooth superlubric sliding with reduced kinetic friction ensues: in such conditions, solid lubrication can provide quite low friction. In practice however neither the substrates nor the lubricant are likely to maintain undefected crystalline order. Defects and/or impurities between the sliding surfaces, even if diluted to a weak concentration, may suffice to induce pinning and finite static friction [195], thereby eliminating superlubricity. In the opposite condition of prevailing lubricant-substrate interactions, the surfaces are covered and protected from wear by lubricant monolayers: sliding occurs inside the lubricant bulk. In such condition the lubricant film can be led to melt during sliding; alternatively, the layering imposed by the surfaces can remain solid, with slips occurring in a layer-over-layer sliding [144,148].

Simulation of extreme frictional regimes.
Simulations are of particular value in the exploration of extreme frictional regimes, that are difficult to access experimentally. Among such extreme regimes, researchers have investigated or are investigating high temperature, high speed, high pressure, and high plate charging in ionic liquid lubrication. Although for most of these conditions there still is no experimental evidence to discuss, simulation has made some interesting predictions that should become of future reference.
High temperature. Close to the substrate melting point T m , its crystal surface may or may not undergo surface melting -the formation, in full thermal equilibrium, of a thin liquid or quasi-liquid film at the substrate-vacuum interface [196]. Either event influences importantly the contact of an AFM tip with the surface. Surface melting gives rise to a local jump-to-contact of the film with the AFM tip, as was found both in experiments [197] and in MD simulations [198]. In that case, friction is expected to become hydrodynamic and uninteresting. More interesting is the case where the substrate surface does not undergo surface melting, such as is the case for particularly well packed, stable surfaces like Pb (111) or NaCl(100) [196]. For an AFM tip sliding on NaCl(100), frictional MD simulations suggested two quite different outcomes depending on the frictional mode. A sharp penetrating tip plowing the solid surface experiences a large friction, which drops sharply when the substrate temperature is only slightly below T m , so that the Joule heat suffices to raise temperature locally and form a liquid drop accompanying and lubricating the moving tip. A blunt tip sliding wearlessly experiences instead a very small friction at low temperature, counterintuitively surging and becoming large close to T m , where the nonmelting surface lattice softens -a phenomenon analogous to that exhibited by flux lattices in type II superconductors [199].
High speed. Friction at high speed, of order of 1 m/s or higher, is common in several technologically relevant situations but is rarely addressed in nanoscale, atomistically characterized situations, where velocity is more typically 1 µm/s, many orders of magnitude smaller. As anticipated in Sec. 4.2, MD simulation is an ideal tool for the study of friction in fast-sliding of nanosized systems. Using gold clusters on graphite as test system, simulation has explored high-speed friction, and especially differences and similarities from low speed, examining the slowing down of a ballistically kicked cluster. Both kinetic frictions are similarly viscous -proportional to velocity. However, they show just the opposite thermal dependence. Whereas low speed (diffusive) friction decreases upon heating, when diffusion increases, the high speed (ballistic) friction rises with temperature, when thermal fluctuations of the contact increase [200].
High pressure. The local uniaxial pressure transmitted to a local contact by the overall load on a slider may reach a hundred kbar, but is generally not very well characterized, and the effects of pressure insufficiently explored. MD simulation makes suggestions of different kinds. First, pressure may provoke structural transformation of a crystalline substrate (or slider) from its initial crystal structure to another. As a recent simulation has shown [112] this will reflect in a frictional jump, either up or down. Second, pressure may bring a solid compound close enough to the chemical stability limit for the frictional perturbation to cause bond breaking and the beginning of chemical decomposition [201]. Third, pressure may lead to electronic or magnetic transformations, such as insulator-metal transitions, and this may also in principle influence friction.
High plate charging in ionic-liquid lubrication. Ionic liquids are salts whose ions have such a large size that T m falls below room temperature. Experimental data have shown that friction across contacts lubricated by ionic liquids depends on the state of electrical charging of the sliders [202,203]. MD simulations applied to heavily simplified ionic liquid models indicate how this dependence can be ascribed to electrically induced structural modifications at the slider-lubricant interface [204][205][206]. For extreme plate charging, these modifications may even modify the lubricant thickness, and also affect its whole molecular structure, with strong predicted consequences on friction [207].

Earthquakelike Models
On meso-and macro-scales the interface between two bodies is quite generally far from uniform. When rough surfaces come into dry contact, the actual contacts occur at asperities of different sizes, typically characterized by a fractal distribution. Even for a contact of ideally flat surfaces of polycrystal bodies, different regions of the interface will be characterized by different local values of the static friction due to structural or orientational domains [111,191]. For a lubricated contact, different values of the local static friction may appear due to patches of solidified lubricant or due to capillary bridges. All these cases can be rationalized with the help of an earthquakelike (EQ) model first introduced by Burridge and Knopoff [209] to describe real earthquakes. The nature of two problems -earthquakes and friction -is very similar: differences are restricted to their spatial-temporal scales: kilometers and years to millenia in geology compared to nanometers and seconds to hours in tribology. The EQ model, also known as the spring-and-block model or the multi-contact model, has been successfully used in many studies of friction [141,[210][211][212][213][214][215][216][217]; similar schemes have been used also to model the failure of fiber bundles and faults [218][219][220].
In EQ models, two corrugated surfaces make contact only at a discrete set of points, as shown schematically in Fig. 13. When the slider moves, a single point contact elongates elastically as a spring, as long as the local shear force f i = kx i (x i is the contact stretching and k is its elastic constant) remains below a threshold value f si = kx s ; then the contact breaks and slips for some distance, as indeed was observed in tip-based microscopy experiments as well as in MD simulations. Then, either immediately or after some delay time, the contact is reformed again with zero or lower stretching and a new threshold value.
The simplified version of the EQ model assumes that all contacts have the same threshold f si ; such a model however corresponds to a singular case and may lead to unphysical results [221]. In a real situation, the contacts always have different thresholds with a continuous distribution P c (x s ) of their static threshold elongations. Therefore, when the upper block begins to advance, the forces acting locally on each contact increase: at successive moments, the contacts begin to snap in a sequence: weaker contacts break earlier, while the strongest contact resists to the last.
EQ-like models are usually studied by simulation. Nonetheless, the kinetics of the EQ model can be described by a master equation (ME), occasionally known as Boltzmann equation or kinetic equation. In concrete, indicate the distribution of contact's stretchings x by Q(x; X), when the sliding plate reaches position X (see Fig. 13): the evolution of Q(x; X) is described by the equation [208,221] where P (x)∆X is the fraction of contacts which break at stretching x as a consequence of the plate advancing by ∆X. The "rate" P (x) and the distribution of the breaking thresholds P c (x) are connected by the relation indicating that the fraction of contacts which snap when X increases by ∆X equals those that have their thresholds between x and x + ∆X, divided by the total fraction J c (x) of contacts still unbroken at stretching x.
The EQ model can be extended to account for thermal effects as well as the ageing of contacts [221]; the latter requires an additional equation to describe the increase of threshold values with the time of stationary contact. Analytic solutions of the ME are available and, in the smooth-sliding regime, they provide us with the velocity and temperature dependence of the kinetic friction force F k [222]. Contrary to the Amontons-Coulomb laws, which state that (macroscopic) friction is independent of velocity, the friction force in the EQ model depends on the sliding speed. At small driving velocity F k increases linearly with speed, F k (v) ≈ γ * v. Indeed, if the slider moves slowly enough, thermal fluctuations will soon or later break all the contacts. The slower the slider the longer time any contact will be given to undergo a fluctuation exceeding its respective threshold; therefore the smaller the resulting kinetic friction force is. This linear F k (v) dependence could be represented as a (characteristically large) effective viscosity γ * of a ultrathin lubricant film [214,216]. At high driving velocities the friction force exhibits the opposite behavior, it decreases when v grows, F k (v) ∼ f s v * /v, due to relaxation, which one can also call an ageing effect. After snapping, a contact slips for a short time, then it stops and is reformed again, growing in size and strength. The faster the slider moves, the shorter time the contacts are left to be reformed and to grow. Overall therefore, the kinetic friction force F k (v) increases with v at low v, up to a maximum at v 0 ∼ (f s v * /γ * ) 1/2 and then decreases. At very high velocities F k (v) should eventually grow again due to an increased damping in the slider bulk. At intermediate speeds typical of experiments, the interplay of thermal and ageing effects generates a weak (approximately logarithmic) F k (v) dependence, approximately consistent with the Amontons-Coulomb laws, although the proper F k (v) dependence is rather difficult to detect in experiments [210,211].
On the decreasing ("velocity-weakening") branch of the F k (v) dependence at v > v 0 the slider motion may become unstable and change from smooth sliding to stickslip motion. If the slider velocity increases due to a fluctuation, the friction force decreases, and the slider accelerates. This effect is studied usually with the help of phenomenological approach (e.g., see [144] and references therein). A detailed study of the EQ model [214,216,221,223] shows that stick-slip for the multi-contact interface may appear, if and only if two necessary conditions are satisfied. First, the interface must exhibit an elastic instability. When the slider moves, the contacts break but then are formed again, but only if the reformed contacts build up a force sufficient to balance the driving force, the motion will be stable. Otherwise the slider will develop an elastic instability, and will keep accelerating until the overall pulling spring force (of elasticity K, see Fig. 13) decreases enough to regain stability. Second, the contacts must undergo ageing. Once these conditions are satisfied, stick-slip will exists for an interval of driving velocities only, v 1 < v < v 2 , while for lower (v < v 1 ) and higher (v > v 2 ) speeds the motion is smooth [223].
The ME approach discussed above however assumed a rigid slider which is not a proper model of a realistic extended system. For a nonrigid slider, its elasticity produces contact-contact interaction: as soon as a contact fails, the forces on nearby contacts must increase by an amount δf . This δf (r) was shown [224] to depend on the distance r from the failed contact as δf (r) ∝ r −1 at short distance r ≪ λ, and as δf (r) ∝ r −3 at long distance. The crossover length λ ∼ a 2 E/k, known as the elastic correlation length [224,225], depends on properties of both the slider (its Young modulus E) and the interface (the mean separation a between nearby contacts and their average rigidity k). Accordingly, a simpler model can be formulated that considers the slider as rigid across regions of lateral size ∼ λ, with the micro-contacts inside each λ-sized area treated as a single effective macro-contact. For this "λ-contact" the parameters can be evaluated by solving a specific ME.
Numerics also showed, in addition, that a large fraction of the extra intercontact force concentrates behind and in front of the snapped contact, implying that effectively the interface can be treated as a one-dimensional chain of λ-contacts, at least approximately. Now, if a λ-contact undergoes elastic instability, namely if at a certain threshold stress the λ-contact snaps and then advances, then the surrounding λ-contacts acquire an extra chance to also fail: this mechanism results in a sequence of snaps propagating forward and backward along the interface as in a domino effect. The resulting dynamics of the chain of λ-contacts could then be described by the Frenkel-Kontorova model of Sec. 3 above, but replacing the sinusoidal substrate potential with a sawtooth profile, i.e. a periodically-repeated array of inclined lines [226]. This approach allows one to calculate the maximum and minimum shear stress for the propagation of this self-healing crack (the minimum shear stress coincides with the Griffith threshold [227]), and also the dependence of the crack velocity on the applied stress. When a λsized contact fails at some point along the chain subject to a uniform shear stress, two self-healing cracks leave the initial snapping contact, propagating in opposite directions as divergent solitons similar to the kink-antikink pair of Fig. 2, until these cracks either reach the boundary or meet another crack generated elsewhere.
On the other hand, nonuniform shear stress is relevant for experiments such as those carried out by Fineberg's group [228][229][230] with the slider pushed at its trailing edge: at this location the shear stress is maximum and across the block it decreases with the distance from the trailing edge. In this system, the leftmost λ-contact is the most likely to fail first, as the pushing force is increased. This failure will result in an increased stress concentrating on the successive λ-contact, which will fail as well. This process will repeat itself until a self-consistent stress remains below the breaking threshold everywhere in the slider. As a result, the self-healing crack initiated at the trailing edge will propagate along the interface for a certain distance Λ (that can be calculated [231]), releasing the stress at its tail side, while accumulating extra stress at its forward side. When the pushing force further increases, a second crack starts at the trailing edge, and can trigger a failure sequence in the pre-formed stressed state, thus propagating to some extra distance. These multiple cracks repeat themselves until they reach the slider leading edge, resulting in a major collective slip. Thus, at the sliding onset, several cracks advance along the interface, with the whole slider undergoing multiple small forward slips (the so-called precursors). In experiment these precursors could be detected and could help predicting the eventual large "earthquake" [215,230].
As mentioned at the beginning of this section, the EQ model was formulated initially to explain earthquakes. Actual earthquakes follow two approximate empirical laws: the one named after Gutenberg-Richter [232,233] states that the number of earthquakes with the magnitude ≥ M scales with the corresponding magnitudes according to a power law; the Omori law [234] states that aftershocks occur with a frequency decreasing roughly with the inverse time after the main shock. EQ-like models discussed above can provide a rationale for both these laws of seismology. Specifically, the Gutenberg-Richter law can be understood as a direct consequence of contact ageing [235]; the Omori law can be interpreted in terms of cracks propagating to a finite distance -after a major earthquake, the stress is not released in full: a certain amount of stress remains stored at a distance ∼ Λ from the main shock [154], where an aftershock is likely to occur some time later.

Conclusions
The fascinating and multidisciplinary topic of microscopic friction, where physics, engineering, chemistry and materials science meet to study the process of converting mechanical energy irreversibly into heat, still lacks fundamental understanding, and increasingly calls for well-designed experiments and simulations carried out at wellcharacterized interfaces. Although AFM, SFA, and QCM setups are providing insight in the high nonlinear out-of-equilibrium interface processes at the small length scales, these advanced experimental techniques still provide averaged tribological data. Overall physical quantities, such as the average static and kinetic friction, the mean velocity and slip times, do not allow to tackle easily the problem of relating the mesoscopic frictional response of a driven system to the detailed microscopic dynamics and structural rearrangements occurring at the confined interface under shear.
In this respect, by explicitly following and analyzing the dynamics of all degrees of freedom at play in controlled numerical "experiments" with interface geometry, sliding conditions, and interparticle interactions can be tuned, mathematical modeling and computer simulations have proven remarkably useful in the investigation of tribologic processes at the atomic scale and are likely to extend their role in future frictional studies.
Even though a large number of open questions remains to be addressed, these modeling frameworks have provided effective insight in the nonlinear microscopic mechanisms underlying the complex nature of static and kinetic friction, the role of metastability, of crystalline incommensurability, and of the interface geometry. Each theoretical approach, from simplified descriptions, to extended realistic MD and hybrid multiscale simulations, has limitations and strengths, with specific abilities to address specific aspects of the physical problem under consideration. Thus, a robust prior understanding of the theoretical background is a basic first step in deciding which modeling features deserve a specific attention and which ones are rather irrelevant details, and then in selecting the best methodological approach for a given problem.
Concluding, it is worth recalling that novel experimental approaches address the intrinsic tribological difficulties of dealing with a buried interface with a very limited control of the physical parameters of the frictional system: artificial systems consisting in optically trapped charged particles, either cold ions in empty space or colloidal particles in a fluid solvent, forced to slide over a laser-generated periodic potential profile. Indeed, especially in 2D colloid sliding it is possible to follow each particle in real time, like in MD simulations. By knowing and, on top of that, tuning the properties of a sliding interface, our physical understanding can expand significantly and open up possibilities to control friction in nano-and micro-sized systems and devices, with serious possibilities of bridging between nanoscale and mesoscale sizes and phenomena.