The breakdown of superlubricity by driving-induced commensurate dislocations

In the framework of a Frenkel-Kontorova-like model, we address the robustness of the superlubricity phenomenon in an edge-driven system at large scales, highlighting the dynamical mechanisms leading to its failure due to the slider elasticity. The results of the numerical simulations perfectly match the length critical size derived from a parameter-free analytical model. By considering different driving and commensurability interface configurations, we explore the distinctive nature of the transition from superlubric to high-friction sliding states which occurs above the critical size, discovering the occurrence of previously undetected multiple dissipative jumps in the friction force as a function of the slider length. These driving-induced commensurate dislocations in the slider are then characterized in relation to their spatial localization and width, depending on the system parameters. Setting the ground to scale superlubricity up, this investigation provides a novel perspective on friction and nanomanipulation experiments and can serve as a theoretical basis for designing high-tech devices with specific superlow frictional features.

In the emerging field of nanoscale science and technology, understanding both the statics and the non-equilibrium dynamics of systems with many degrees of freedom in interaction with some external potential, as is commonly the case in solid state physics, is becoming a crucial issue. Friction belongs to this category too, because the microscopic asperities of the mating surfaces may interlock and give rise to intriguing length scale competition mechanisms. Frequently, despite the basic level of details, simple phenomenological models of friction 1,2 have revealed the ability of describing the main features of the complex microscopic dynamics, providing good qualitative agreement with experimental results on nanoscale tribology, or with more complex simulation data of sliding phenomena. With respect to this, the application of driven Frenkel-Kontorova-like (FK) systems 3,4 , modeling the dissipative dynamics of a chain of interacting particles that slide over a rigid substrate potential due to the application of an external driving force, has found increasing interest as a possible interpretative key of the atomistic processes occurring at the interface of two materials in relative motion.
In particular, the remarkable idea of frictionless sliding connected with interface incommensurability, a pervasive concept of modern tribology named superlubricity 5 , can be mathematically drawn in the framework of the FK model 6 . When two contacting crystalline workpieces are incommensurate, due to lattice mismatch or to angular misalignment, the minimal force required to achieve sliding, i.e. the static friction, should vanish, provided the two substrates are stiff enough. In such a geometrical configuration, the surface mismatch can prevent asperity interlocking and collective stick-slip motion of the interface atoms, with a consequent negligibly small frictional force.
Systems achieving low values of dry sliding friction are of great physical and potentially technological interest, e.g., to significantly reduce dissipation and wear in mechanical devices functioning at various scales. Superlubricity, with usually consequent ultra-low dynamic friction, is experimentally rare, and has been demonstrated or implied in a relatively small number of cases, such as sliding graphite flakes on a graphitic substrate [7][8][9] , cluster nanomanipulation 10,11 , sliding colloidal layers [12][13][14] , and inertially driven rare gas adsorbates [15][16][17] .
Until very recently, superlubricity has been practically observed only on the nanoscale. A short time ago, a breakthrough has been achieved, demonstrating the existence of superlubric regime of motion for micrometer size graphite samples 18 and centimeter-long double walled carbon nanotubes (DWCNTs) 19 . With the current capability of synthesizing and manipulating quasi 1D atomically perfect objects of extended length, such as telescopic nanotubes 19 , graphene nanoribbons 20 , aromatic polymers 21 , or soft biological filaments 22 , nanotechnologies open now the possibility to transpose the peculiar nanoscale tribological properties to larger scales and exploit them to control sliding-induced energy dissipation in state-of-the-art technological devices.
Though, the robustness of the superlubricity phenomenon remains a challenge, and the conditions of its persistence or the mechanisms of its failure are cast as key questions to be addressed. Even in clean wearless friction experiments with perfect atomic structures, superlubricity at large scales may surrender due to elasticity of contacting samples 23 .
A recent letter 24 , via a FK modeling approach, has shown that the effective rigidity of an incommensurate slider, necessary condition for the occurrence of a superlubric regime, is strongly affected by the driving protocol used to induce motion. In an edge-driven configuration, where the external pulling or pushing force is applied to the edge of the slider, and which is typical for many tribological and nanomanipulation experiments performed with atomic force microscope, a slider critical length (depending on intrinsic material properties and experimental parameters) emerges above which superlubricity breaks down. Here, we highlight new features of the transition from superlubric to high-friction state which occurs above the critical length. In particular, we consider different driving and commensurability configurations, discovering the occurrence of previously undetected multiple jumps in the friction force as a function of the slider length. We analyze in details the mechanisms leading to the striking boost in dissipated energy materializing abruptly, above the critical size, due to a localized deformation of the slider becoming commensurate to the underneath substrate lattice and thus giving rise to dynamical high-friction stick-slip regimes.
The paper is organized as follows. In the next section, we present the main features of the edge-driven FK-type model, briefly summarizing the simulation procedure. The two subsequent sections are devoted to the numerical results elucidating the novel mechanisms of energy dissipation leading to the breakdown of superlubricity for long extended systems, together with a consistent analytical framework for the determination of the length critical size. Section on domain wall features discusses relevant characteristics of the driving-induced commensurate dislocations in the slider, specifically in relation to their spatial localization and width depending on the system parameters. The novel occurrence of multiple commensurate walls is also presented. The penultimate section deals with the evaluation of the critical length hampering superlubricity in practical experimental systems, such as DWCNTs. Conclusions are given in the last section.

The model
The mechanical and tribological properties of the quasi 1D nanostructures previously introduced can be analyzed in the framework of the Frenkel-Kotorova model 3,4 . Specifically, they are treated as finite chains of N particles of mass m, connected by springs of stiffness K, having rest length a c , and driven on a sinusoidal potential with periodicity a s and amplitude U 0 representing the interaction with the substrate. As sketched in Fig. 1(a) the rightmost particle of the chain, with coordinate X N (t), is pulled or pushed at constant velocity V 0 through a spring of constant K dr representing the lateral stiffness of the AFM cantilever. As in real AFM experiments the friction force F(t) is estimated by the lateral deflection of the cantilever, i.e. by the spring elongation F = k dr [V 0 t − X N (t)]. A useful dimensionless unit system is obtained setting a c as the unit length, τ = / m K as the unit time and Ka c 2 as the unit energy, with this choice the equations of motion are: with the substrate force given by In the dimensionless unit system 0 (a c /τ) and k dr = K dr /K. The work done by the external force excites the internal degrees of freedom of the driven nanostructure first, e.g. molecular modes, and is subsequently dissipated through the substrate degrees of freedom, i.e. phonons. To dispose off the extra energy injected by the external driving a viscous damping term has been included in the equations, γ = ητ is the dimensionless viscous coefficient if η is the dimensional one, 1/η represents the characteristic time necessary for the energy dissipation.
The equations of motion (1) have been integrated numerically using a velocity-Verlet algorithm. Finite temperature simulations have been also performed by means of a Langevin thermostat, however we found that up to T = 0.4U 0 /k B the thermal effects do not lead to qualitative changes, thus we focus here on the case of T = 0 only.
All the quasi 1D systems previously discussed present some common feature: (i) their periodic atomic structure is incommensurate with respect to the substrate one; (ii) their interatomic stiffness K is sufficiently larger than the interfacial stiffness K int = U 0 (2π/a s ) 2 . These conditions, in the thermodynamic limit (N → ∞), give rise to a continuum set of ground states that can be reached adiabatically through nonrigid displacement of chain atoms at no energy cost (Goldstone mode), with a consequent vanishing static friction 6 .

Breakdown of Superlubricity
Before discussing in detail the novel mechanism of dissipation leading to the breakdown of superlubricity for long chains, we start analyzing the frictional response as a function of the chain length first consid-  In this super-low friction regime the friction force per particle slowly decreases with N, as illustrated in Fig. 1(c). At a critical particle number  N 750 cr a sudden increase of the average friction force breaks the super-low friction regime. The time evolution of the friction force for a chain with N = 1000 > N cr is reported in Fig. 1(b) showing an abrupt jump during the initial transient almost doubling the average friction value in the steady sliding state plateau. A quantity intimately connected to friction is the energy dissipation, to better understand which changes affect the sliding motion for N > N cr one can analyze how the power dissipated by the viscous damping is distributed along the chain. More precisely, in our dimensionless unit system, one can plot a normalized time-averaged power dissipated by the i-th atom: where, in the time interval T over which the average is taken, the chains already reached the steady sliding state. For N < N cr such a quantity is almost zero everywhere except at the edges, due to a slightly higher particle mobility but, as shown in Fig. 2(a), for N > N cr a narrow region appears, in a specific point of the chain, where dissipation increases by two orders of magnitude. A different chain arrangement must evidently give rise locally to a completely different and extremely dissipative sliding mechanism.
To further elucidate the local rearrangement of the particles it is convenient to plot the local commensuration index: i.e. a quantity close to one when two neighboring particles are locally in registry with the substrate and larger or smaller than one if the two particles are closer or far apart than a s . A color map representing the time evolution of the commensuration index, from the initiation of driving to the onset of the steady . sliding state, is plotted in Fig. 2(e). A narrow yellow region where nearby particles are commensurate ( ( )  d i 1) nucleates at the driving edge and propagates backward in the chain upon reaching a certain constant position L cr . It works as a domain wall separating two incommensurate chain portions with d(i) < 1 on the left and d(i) > 1 on the right. Figure 2(i) shows a cross section of the 2D color map taken at a large simulation time where the steady state is already reached, it clearly shows that the chain undergoes a linear stretching with a discrete jump in correspondence of the domain wall. It is easily understandable, and readily verified looking directly at the particle trajectories, that stick-slip motion insets within the wall region with a consequent enhancement of dissipation and loss of superlubricity.
Changing the driving mode from pulling to pushing, a jump in the friction force is also observed. As shown in Fig. 2(f) this corresponds to the nucleation, at the trailing edge, of a narrow domain wall region in which ( ) .  d i 0 5, i.e. a region of local commensuration where the compressed particles assume a periodicity half of the substrate lattice spacing a s . As in the pulling case this domain wall region, separating two incommensurate domains having d(i) < 0.5 on the left and d(i) > 0.5, propagates forward in the chain up to a certain location L cr . Figure 2(b) shows the occurrence of a dissipation peak in analogy to the pulling case. The overall friction force experienced by the chain is smaller than in the pulling counterpart, while a 1/1 commensuration requires in fact the wall particles to sit in the very bottom of the surface potential, with a consequent large effort to make them sliding forward, in a 0.5/1 commensuration the two particles sharing the same potential valley will lie far from the minimum, less force is needed to push the wall forward.
Up to now we only considered overdense chains but, given the asymmetry between pushing and pulling driving, it is definitely interesting to investigate also the behavior of underdense chains, i.e. the case . Also in this cases pushing and pulling give rise to very different behaviors as shown by Fig. 2(g,h). Again we have the nucleation and propagation of domain walls where the particles get commensurate with respect to the substrate potential, however, given the lower density of particles in the chain, upon pulling we reach a 2/1 commensuration while upon pushing the particles are squeezed back to a 1/1 configuration, see Fig. 2(k,l). In both cases the commensurate particles live in the very bottom of the surface potential and, as in the overdense pulling case, the total friction force jumps by more than a factor two for N > N cr . Looking at the average dissipated power, shown in Fig. 2(c,d), one might notice again a large peak in correspondence of the domain wall region, but also a pronounced dissipation elsewhere, contrary to the overdense analogs. In the underdense case this can be explained by a much larger variety of partially commensurate slider-substrate configurations, leading to a broader multiplicity of dissipation peaks.

Analytic Treatment
In this section we derive an analytic expression for the critical number of particles in the chain N cr above which superlubricity is broken. The explicit calculation will be carried out for the pulling of overdense chains only, the same arguments can be straightforwardly applied to the other three cases above presented.
We start writing the time-average friction force for a generic chain in its steady sliding state 25 : This equation is exact and is derived simply imposing energy conservation, i.e. the energy pumped into the chain by the external driving per unit time equals the energy dissipated per unit time by the viscous damping applied to the particles For chains with N < N cr , the absence of stick-slip allows us to apply a linear perturbation theory writing is a small dimensionless parameters representing the ratio between the chain and substrate stiffness, =  x v i 0 0 is the exact solution in absence of the substrate potential, i.e.  = 0. Substituting the velocity expansion in eq. (6) we get 25 all the unknowns hidden inside α 1 , while it remains hard to get a simple analytic expression for this parameter, we can easily find its numerical value by fitting the simulation results. To this aim we performed many different calculations of the average friction force varying U 0 , K, V 0 and η for chains with number of particles N b < N < N cr where 〈 F〉 /N is almost independent of N, i.e. in the initial plateau of Fig. 1(c). Choosing α 1 = 2.64 ± 0.11 we obtain a remarkable agreement between the analytic expression (8) and the numerical results with  . < < . 0 03 03 as shown in Fig. 3(a). Only for large U 0 , namely  = . 0 75, we start to see a significant discrepancy between the model and the simulations, for such a large  value the adopted linear expansion is clearly insufficient and higher order therms should be included in eq. (8).
Once again we stress that eq. (8) describes a viscous-like friction, proportional to V 0 and η and only weakly dependent on U 0 and K and is suitable for N < N cr . For N > N cr the domain wall region nucleating inside the driven chain is characterized by a stick-slip type of friction, increasing with U 0 and K and practically independent of V 0 and η, such an highly non-linear kind of motion cannot be approximated with a power expansion. Our simulations show that above the critical length the friction force can be calculated as a sum of the contributions given by eq. (8) and by the dissipation in the local region at the domain wall.
As discussed in the previous section the maximum length at which the driven chains exhibit a super-low friction behavior is limited by the nucleation of a localized commensurate region at the driving edge. The external force F ext necessary to commensurate the N and N − 1 particles with the underneath substrate potential can be easily calculated obtained when the N-th particle sits on a maximum. Clearly for the other three cases presented in Fig. 2 the nucleation condition, to be imposed in order to calculate F ext cr , is different as the domain wall region is constituted by a different arrangement of the particles, as schematized by the cartoon insets of panels (a)-(d). Expression (10) is exact and can be directly verified through the simulations. As shown in Fig. 4 the F ext cr value calculated from eq. (10) for different K and U 0 , dashed lines, always match the value at which F(t) jumps vertically in the corresponding simulations, i.e. the force value at which the commensurate region nucleates. Moreover, according to eq. (10), the critical force does not depend on the length of the chain and on the dynamical parameters, i.e. the viscous damping and the driving velocity. This is also easily verified through our simulations, Fig. 5 shows again F(t) during the nucleation process for different η, V 0 and N and is clearly seen that, although the nucleation takes place at different times, F ext cr is constant. Now, in a steady sliding regime the time-averaged pulling force equals the average friction force and an expression for the critical number of particle N cr can thus be obtained equating (8) and (10) ]. This new expression shows that the critical number of particles increases with the stiffness of the chain and decreases with increasing the damping coefficient and pulling velocity. N cr also decreases with increasing the ratio between the two stiffness, K int /K, however this effect is weak since K int /K is a small parameter. To verify the robustness of our model we tested the prediction of eq. (11) against the numerical simulations performed with a broad range of parameters V 0 , η, U 0 and K, the results are shown in Fig. 3(b). For all values of system parameters the theoretical results agree well with numerical simulations. This result show that the critical number of particles is mostly determined by the stiffness of the chain K, the damping coefficient η, and pulling velocity V 0 rather than by the ratio K int /K as suggested in previous theoretical works 23 . In this work, to magnify and highlight the breaking of superlubricity, we only considered the incommensurability of the chain and substrate structures corresponding to the golden ratio, for which the difference between two periods is relatively large. In real systems with such a large period mismatch anharmonic effects are expected to play a role, but we believe they do not alter the overall picture obtained within the linear elasticity approximation adopted in our numerical and analytic models. Moreover, the results here discussed are valid also for smaller misfits between the contacting lattices. Then the stretching of the chain at the transition is significantly smaller and anharmonic effects can be neglected again. Finally it should be noted that our calculations of N cr are based on the assumption that all particles of the chain experience the same microscopic friction proportional to the viscous damping coefficient. However, in realistic systems there could be additional contributions to friction coming from chemical interactions between the edge particles and the surface and/or defects in the chain structure. These effects will lead to an increase of the average friction force thus reducing the critical number of particles, which could benefit the experimental investigations of the mechanisms discussed here.
Repeating the same fitting procedure for the underdense condition of / = ( + )/ a a 1 5 2 c s , we find α 1 = 2.69 ± 0.15 the same value as before within the estimated error. Figure 3(c,d) report the comparison

Domain Wall Features
As discussed in section 1 and immediately visible in Fig. 2, upon reaching the steady sliding state, the narrow domain wall region takes a well defined and constant position L cr along the chain. This position is determined by the chain and substrate stiffness and by the dynamic parameters m, η, v 0 . On the contrary L cr is not affected by an increase of the chain length, upon increasing N the distance between the free edge and the domain wall is unchanged while the distance between the domain wall and the driven edge increase linearly. The only other quantity displaying such a behavior is the chain stretching, i.e. the slope of the curves in Fig. 2(i-l), influenced by substrate and chain stiffness and by the driving conditions but independent of N both above and below N cr . It is thus clear that, once the driving dependent stretched state is reached, L cr is solely determined by the interplay between the periodic substrate potential and the linearly increasing periodicity of the chain. The position of the i-th particle in a linearly stretched chain is given by X i = a c (i − 1) + δ(i − 1)i/2, where the slope of linear stretching δ = δ(k, v 0 , η, U 0 , m) sets the overall length of the stressed chain, once the steady sliding state in reached, as a function of the driving parameters. The potential and force felt by every particle due to the interaction with the substrate potential is thus readily obtained:  Fig. 6(a) for increasing δ values in the overdense puling case: upon increasing the stretching of the chain a plateau region appears, i.e. a region where nearby particles feel the same substrate potential meaning that they are almost perfectly commensurate. Keeping increasing δ the commensurate region moves backward, in our simulations this backward motion stops when the steady sliding state is reached and no further stretching/elongation of the chain occurs. Elsewhere in the chain the potential felt is rapidly and incoherently oscillating, upon averaging over few neighboring particles the resulting force will thus be zero everywhere except around the commensurate region as show in Fig. 6(b). An expression for N cr can thus be obtained as a function of δ for the pulled overdense chain imposing X N − X N−1 = a s , this yields: From the simulated chain length in the stretched steady sliding state we can estimate δ and use it the two previous equations to evaluate N cr or L cr . The red curve in Fig. 6(a) has been obtained with δ = 8.35 × 10 −4 corresponding to the chain of Fig. 2(a), indeed the position of the plateau in the potential coincides with the position of the domain wall, the estimated critical number of particles N cr = 741 is in good agreement with the simulations, see Fig. 1(c). One might be tempted to compare eqs. (14) and (11) to obtain an analytical expression for δ(k, v 0 , η, U 0 , m). Figure 7 shows the comparison between such a theoretical prediction and the numerical values coming from the simulations. The agreement is fairly good considering that eq. (14) is exact while eq. (11) comes from a power expansion for small U 0 only. As previously discussed and shown in Fig. 2, the dissipated power P(i) exhibits a high peak localized in the commensurate region of the chain, where it takes a value two orders of magnitude higher than in the rest of the chain. The appearance of this peak is the finger print of transition from the super-low to high friction when N > N cr . Figure 8(a) shows that, in the range of parameters here considered, the width of the dissipation peak scales as / K K int . The particles located in the narrow commensurate region corresponding to the domain wall perform stick-slip motion, and their dissipated energy is proportional to the amplitude of the particle-substrate interaction U 0 . As a result, the integral power dissipated by the particles located in the commensurate region scales as K K int . Figure 8(b) shows that this prediction is again in good agreement with the numerical simulations from our simulations, notice that this integral dissipated power is independent of the viscous damping coefficient η and driving velocity V 0 . The power dissipated by the particles located outside the domain wall is proportional to ηV 0 and only slightly dependent on the parameters K and K int , see eq. (8). Thus, the mechanisms of energy dissipation in the commensurate and incommensurate domains of the chain are completely different, originating contributions to the net friction force that can differ by several orders of magnitude.
If, as just demonstrated, the occurrence of a commensuration region in the chain is due to a local matching between the substrate periodicity and the linearly increasing chain periodicity, it becomes essential to check whether this commensuration occurs only once or periodically along the chain. Figure 6(c) shows the potential energy per particle as δ is increased further after the occurrence of the first plateau, indeed many other plateaus appear and condensate within the chain. One can also fix δ increasing the chain length as in Fig. 6(d) (the average force is plotted in place of the potential in order  to better appreciate small variations on a large N interval) force peaks (i.e. plateaus in the potential) appear regularly as N grows in a sort of Moire pattern with δ dependent periodicity. As a consequence, for chains long enough and under the proper driving conditions, more domain wall regions are expected to nucleate leading to higher discrete jumps in the friction force. This behavior is illustrated in Fig. 9 obtained simulating a chain under the same driving conditions and parameters of Fig. 2(a) but four times longer. Two identical domain wall regions propagate in the chain causing a considerable enhancement of friction and giving rise to localized dissipation peaks due to stick-slip sliding. Notice that the second jump occurs when the value of the pulling forces reaches the critical value The corresponding average force exerted by the substrate potential is given by the red plot in Fig. 6(d), also in this case the position of the average force peaks (read potential plateaus) corresponds nicely with the domain wall location along the chain. A similar treatment could be done for the other three cases discussed in Fig. 2 of pushed and underdense chains, naturally the condition to be imposed to derive eqs. (15) and (14) varies from case to case.

Critical Length Estimation for Experimental Systems
Before concluding we give a practical example of how eq. (11) can be used to estimate the critical length for different systems whose tribological properties are experimentally accessible. Among other, particularly interesting is the case of DWCNTs where super-low friction has been observed up to the centimeter length-scale 19 . For these systems the intrashell and intershell stiffness can be estimated as K = Ea c and K int = Ga s 23 , respectively, where the Young s modulus E ~ 0.5 − 1 TPa 26 and intershell shear modulus G ~ 1 GPa 27 have been measured experimentally. Experimental 28 and numerical 29 estimations of the damping coefficient give η ~ 1 ps −1 . The mass m, can be calculated from the bulk density of carbon nanotubes as m = ρ V πDha c where ρ V = 2240 kg/m 3 is the bulk density, h = 0.34 nm is the wall thickness 29 and D is the diameter. As an example of incommensurate DWCNTs consider for instance an armchair tube inserted into a zigzag one, in such a case a s − a c is roughly 0.1a C−C , where a C−C = 0.142 nm is the carbon-carbon bond length. For DWCNTs with the diameter of the outer shell D ~ 2.73 nm pulled with speed V 0 ~ 1μm/s through a force probe attached to one its ends 19 , L cr is estimated to be 0.5 m. Our theoretical model is thus partially validated by the experiment where the longest DWCNT manipulated was ~9 mm and is found to display a superlubric behavior. Moreover eq. (11) suggests that carbon based nanostructures, with a large internal stiffness K due to the C − C bond strength and a small K int due to the weak interplane interaction, are best candidates for further scaling up of superlubricity.

Conclusions
Summing up, in the framework of a Frenkel-Kontorova-like model, we have addressed the robustness of the superlubricity phenomenon in an edge-driven system at large scales, highlighting the mechanisms leading to its failure due to the slider elasticity. The results of the numerical simulations perfectly match the length critical size derived from a parameter-free analytical model. By considering different driving and commensurability interface configurations, we have explored the distinctive nature of the transition from superlubric to high-friction sliding states which occurs above the critical size, discovering the occurrence of previously undetected multiple dissipative jumps in the friction force as a function of the slider length. These driving-induced commensurate dislocations in the slider are then characterized in relation to their spatial localization and width depending on the system parameters. This investigation provides a novel perspective on friction and nanomanipulation experiments and can serve as a theoretical basis for designing nanodevices with specific superlow frictional features.