Robustness of structural superlubricity beyond rigid models

Structural superlubricity is a theoretical concept stating that the friction force is absent between two rigid, incommensurate crystalline surfaces. However, elasticity of the contact pairs could modify the lattice registry at interfaces by nucleating local slips, favoring commeasure. The validity of structural superlubricity is thus concerned for large-scale systems where the energy cost of elastic distortion to break the lattice registry is low. In this work, we study the size dependence of superlubricity between single-crystal graphite flakes. Molecular dynamics simulations show that with nucleation and propagation of out-of-plane dislocations and strained solitons at Bernal interfaces, the friction force is reduced by one order of magnitude. Elastic distortion is much weaker for non-Bernal or incommensurate ones, remaining notable only at the ends of contact. Lattice self-organization at small twist angles perturbs the state of structural superlubricity through a reconstructed potential energy surface. Theoretical models are developed to illustrate and predict the interfacial elastoplastic behaviors at length scales beyond those in the simulations. These results validate the rigid assumption for graphitic superlubricity systems at microscale, and reveal the intrinsic channels of mechanical energy dissipation. The understandings lay the ground for the design of structural superlubricity applications.


Introduction
The notion of structural superlubricity has recently drawn notable interest for its practical implications in engineering, holding great promises in applications that are friction-free, dissipation-free, and even wearless [1][2][3][4][5]. Recent studies on structural superlubricity suggest that the interface between graphite with lattice misfit features robust structural superlubricity performance, which can be achieved by rigid twist or elastic straining [6], and be maintained even in extreme environmental conditions including high pressures [7], high or low temperatures [8], in vacuum or ambient air [2,9]. Early efforts were made by exploring the amplitude of friction force or mechanical energy dissipation of a graphite mesa sliding on a graphitic substrate, demonstrating the rotation-angle dependence and six-fold symmetry of the potential energy landscape, where graphite flakes are considered as rigid plates and elastic deformation is neglected [2,10]. Theoretically, this process can be captured in the Prandtl-Tomlinson (PT) model [11], which treats the mesa as a particle and the substrate contact through a fixed potential surface defined by the mesa-substrate interaction. The energy barrier against interlayer sliding decreases as the interface is driven out of lattice registry. In the incommensurate limit, forces on atoms cancel out, resulting in a frictionless contact in the state of structural superlubricity [12][13][14].
However, this picture may not hold if the size of graphitic flakes along the sliding direction L exceeds a critical value c L , and the interfacial shear stiffness (~L) becomes comparable to or lower than the twodimensional (2D) stiffness of graphene (the surface layer of graphite) or graphite [15][16][17][18][19]. Local reorientation of graphene layers is identified in this condition, which improves the lattice registry within finite domains, and dramatically alters the state of superlubricity. This behavior can be described by more realistic models such as the two-chain Frenkel-Kontorova (FK) model [20], and the value of c L depends on the nature of contact that can be quantified through the lattice symmetry, lattice constants, and the interaction strength between the surfaces. Our recent study shows that a series of commensurateincommensurate transitions (localized lattice slips) arise during the sliding processes between graphene bilayers, illustrating the tendency to maintain local lattice registry in the low-energy AB or BA stacking order by sacrificing the energy of elastic deformation [21]. Graphite with its crystalline surfaces cleaved and prepared with incommensurate registry has been demonstrated as an ideal model material for structural superlubricity [22,23]. It was actually concerned in the very early study on structural superlubricity of graphitic systems, supposing that for sufficiently large contacts, superlubricity might break down, as the two lattices are not perfectly rigid, and a network of misfit dislocations should form between the two, the motion of which will dissipate energy [10]. However, later experimental studies did not report notable evidences of this issue, and the effects of elastic deformation in the graphite mesa-substrate systems seem negligible for microscale structural superlubricity. The rigidity assumption is thus expected to apply up to the micrometre-and centimetre-scales [2,[24][25][26]. To clarify the vague robustness of structural superlubricity at extended length scales, we performed molecular dynamics (MD) simulations to explore interlayer sliding behaviors between graphene multilayers. The size effects in terms of length and thickness of the graphite flakes as well as the lattice registry at their interfaces are discussed through analysis of the friction force, lattice deformation, interfacial shear and sliding.

Setups of MD Simulations
The graphite mesa (m)-substrate (s) system is modeled in a ribbon-geometry with periodic boundary conditions (PBCs) enforced in the width direction (along y, 2-3 nm), to save the computational costs (Fig. 1). Graphite flakes with N layers in each component are constructed with an identical length of L = 10 or 200 nm, in a perfect ABAB stacking order, and labeled as m1, m2, , mN, s1, s2, , sN, respectively. Other dimensions are modeled with open boundary conditions. The mesa-substrate interfaces are initialized in Bernal stacking for commensurate contact or with a large misfit angle of     21 8 through rigid twist for "incommensurate" contact (strictly speaking, this is a commensurate non-Bernal interface with a much larger superlattice periodicity, , than the lattice constant of graphene,  1 42 a is the length of C-C bonds in graphene. MD simulations were performed using the largescale atomic/molecular massively parallel simulator (LAMMPS) [27]. The second-generation reactive empirical bond order (REBO) potential [28] is used for the lattice dynamics in graphene layers, and the interlayer interaction is modeled through the registrydependent Kolmogorov-Crespi (K-C) potential implemented with a taper function [29][30][31]. To simulate the sliding experiments, displacement (  ) is applied to two opposite ends (two rows of hexagons for each) of the mesa and substrate, respectively (Fig. 1). The loading velocity is set as t m/s to reduce the dynamical effects. Simulations with lower velocities were also performed to assess the rate dependence, which validates main conclusions of the current work ( Friction force of non-Bernal interfaces with reconstructed moiré structures is analyzed by driving a circular mesa (  1 N and 2) with a radius of  50 R nm to slide over the substrate (  2 N , 120 nm × 120 nm, Bernal stacking with PBC enforced). The initial twist of the contact is   0 7 . Following structural relaxation using the conjugate gradient (CG) algorithm with a force tolerance of Å 1 10 eV   , the moiré pattern self-organizes into a triangular pattern after elastic relaxation, where domains of ABB (from the bottom layer to the top one) shrink into nodes connecting ABA and ABA' domains, displaying a 3 C symmetry (Fig. S3 in the ESM). Displacement is applied to all atoms in the mesa with a constant speed of 10 m/s, and the bottom layer of substrate is fixed. This setup partially restricts structural relaxation at the interfaces, and is chosen only to demonstrate the effects of moiré patterns on the friction stress. The friction stress is calculated by summing over transverse forces on atoms in the mesa and dividing it by the area of contact. This definition measures the kinetic friction or the mechanical energy dissipation during the sliding process, and is consistent with the experimental measurement [25]. On the other hand, the peak values of the transverse forces measure the interfacial shear resistance.

Lattice-registry measures
Lattice slip between neighboring graphene layers is analyzed to measure relative in-plane displacement of carbon atoms during deformation ( Fig. S4 in the ESM) [21]. By correlating the coordinates of carbon atoms located in the top and bottom graphene layers at the contact, the projected position of the i-th atom in the top layer is i x is the current position of the reference-atom triplet (j = 1,2,3), and i j k is the projection factor. The coordinates can then be mapped into an unstrained lattice system of the bottom layer as x . The atomic slip length (SL) is defined in the unit of the length of vector i  x . For interfaces with initially Bernal stacking, SL = 0, 0.5, and 1.0 correspond to the AB, SP, and BA configurations, respectively. For three neighboring layers, SL = 0,1,0 or 1,0,1 corresponds to ABA or BAB, and SL = 1,1,1 or 0,0,0 corresponds to ABA' or BAB'. To distinguish ABB, ABA, and ABA', the registry index (RI) is calculated at the interface ( Fig. S4 in the ESM). Each carbon atom is considered as a circle in the basal plane with a radius of . RI is computed through the overlap area between carbon atoms in different graphene layers [23]. To measure local lattice slips, we define the local registry index (LRI) for the A, B subsites. For Bernal (AB or BA) stacking, neighboring LRIs for the C atoms on the subsites are 1 and 0.

Results
In graphitic superlubricity systems at microscale or www.Springer.com/journal/40544 | Friction beyond, the path of motion is determined by the outlines of graphite flakes [2]. Notably, rotational instability was reported for the favor of cohesive energy at the state of Bernal stacking [32]. We thus studied sliding between graphene or graphite (graphene multilayers in our model) interfaces with initial Bernal stacking, along the armchair and zigzag directions. The nature of sliding friction between crystalline surfaces was reported to feature a notable size effect, where localized lattice slip and deformation are nucleated and develop if the strained solitons can be accommodated [31,33]. Our simulation results show that for L = 10 nm, typical stick-slip phenomena are observed for the graphitic contact, where the whole mesa slides over the substrate (Fig. S5 in the ESM), while for L = 200 nm, spatial patterns of the friction characteristics emerge (Fig. 1). Our following discussion will then be focused on the simulation results for L = 200 nm, and the size effect will be discussed through model analysis.
Tensile stress in mesa/substrate (    M S ) and elongation strain of each layer ( L  ) increase with the sliding displacement (  ) before reaching a plateau (Fig. 2). The amplitude of tensile stress displays abrupt jumps due to the relaxation processes, which are absent in strain. Besides of the linearly-elastic responses, there is a single slip occurring at  cr1 . As the graphene layers are fully stretched, overall sliding is activated after  cr 2 . In contrast to typical interfacial elastoplastic responses to shear in the continuum models [34], the evolution of tensile stress and strain against sliding demonstrates stick-slip features after  cr 2 , which are commonly identified for crystalline interfaces [35]. The interfacial slips (out-of-plane dislocations with Burgers vectors defined by the graphene lattice) or strain solitons (soliton-like structural boundaries between domains with different stacking orders) in the graphene layers continuously develop at this stage ( Fig. 1) [21], and the intralayer elastic deformation becomes significant. The lattice slip at the interfaces and development of strained solitons in the graphene layers characterize the reconstruction of interlayer registry, which is the origin of the modification in the friction stress. These results echo observations in crystal plasticity [36], and signal breakdown of the rigid assumption. The elastic and plastic responses measured in the averaged friction stress  f thus correspond to the static and kinetic friction, respectively. The stick-slip processes result in mechanical energy dissipation from the translational motion into lattice vibration or the phonon bath [37].
The cohesive energy at the graphitic interface maximizes at perfect Bernal stacking and is reduced under the displacement load. Balance between interlayer and intralayer interaction at the atomic scale regulates the interlayer registry and strain distribution in graphene layers. Before the collective motion of out-of-plane dislocations or strained solitons is activated, interfacial shear is localized at the two ends of contact between the mesa and substrate, and the maximum shear increases with the force applied. As the interfacial shear  reaches a critical value  cr , where  L reaches its first critical value  cr1 , local slips or out-of-plane dislocations are nucleated, separating commensurate domains that stay in the Bernal stacking by incommensurate domain walls (Fig. 1). The appearance of domain walls at the interface corresponds to the rise of localized strain in the graphene layers, or namely the creation of strained solitons. As the amplitude of load (end displacement) increases further,  L reaches a second critical value  cr 2 . More strained solitons are then created, growing from the ends of contact and moving towards the centers. These solitons annihilate during collision if the selection rule is manifested [21]. After local slips develop across the whole graphene layers, overall sliding is activated. The friction stress increases with the size of contact L for more local slips are accumulated. The average friction stress  f is measured to be      2 (3 5 0 23) 10 GPa for    1 N L 200 nm and      2 (7 1 1 8) 10 GPa for    5 200 N L nm (Fig. 2). These values are much reduced from the values calculated from rigid models (0.85, 0.83, and 0.18 GPa for sliding along the armchair, zigzag, and AB-SP-BA paths). The thickness dependence of  f in the deformable model is resulted from the increase in the equivalent stiffness of the graphite flakes with the value of N.
Notably, the birth, evolution, and annihilation of strained solitons appear not only for the interface at the mesa-substrate contact, but also well into the graphene multilayers as a result of the load transfer between graphene layers (2-3 layers, see the m1-m2 | https://mc03.manuscriptcentral.com/friction or s1-s2 interfaces for  5 N , Figs. 2(e) and 2(f)). The depth increases with the effective stiffness of graphite, or the number of graphene layers, N. This nonlocality across the graphene layers opens additional channels of mechanical energy dissipation that are overlooked in simplified theories such as those in the PT and two-chain FK models. Beyond the critical strain required for collective soliton development  cr 2 ( ), layers m1 and s1 are much more stretched than the layers m2-mN and s2-sN. The Burgers vectors of dislocations in the mesa are opposite to those in the substrate, demonstrating the centrosymmetry of the simulation setup (Fig. 2(f)). From the spatial distribution of SL (Fig. 3), we find that lattice slips between m1 and s1 for  5 N and the evolution (propagation, annihilation) of strained solitons in the graphene layers behave similarly to that in the graphene bilayer (  1 N ). However, slips between m1 and m2 (s1 and s2) grow from the free ends, and the strained solitons propagate to the other ends without collision or annihilation since the Burgers vectors of the out-of-plane dislocations are in the same direction. Pileup of dislocations at the m1-m2 and s1-s2 interfaces results in the hardening effects (Figs. 2(d) and 2(e)). The stacking order of s2-s1-m1 changes from ABA to BAB' or BAB to ABA' if lattice slips are activated at the s1-m1 interface, and BAB' to ABA' to BAB if slip also occurs at s1-m1 and s1-s2, modifying the load-transfer capability at the interface and the friction stress. Similar phenomenon has been observed for load applied along the zigzag direction, but slips only occur at the contact as the critical strain is much lower (Fig. S6 in the ESM). For sliding along other directions, local AB and AA regions are created, demonstrating the elastic effects.
The length scale under discussion is limited by the capability of MD simulations, which can be extended by analysis using theoretical models. The deformable tension-shear (DTS) model [38] with constant interlayer shear modulus and FK model [21,39] with a trigonometric form of interlayer interaction can capture the development of strained solitons, as well as shear and dislocation dynamics at the interface. Key features including the width, period of strained solitons, and values of  cr1 and  cr 2 can be estimated from the competition between intra-layer elasticity and registrydependent interlayer interaction [21]. From the twochain FK model, the characteristic width of solitons Nevertheless, increasing the thickness of graphite improves the in-plane rigidity and reduces the critical strain. The values of  cr1 and  cr 2 decrease with N ( Fig. 3(a)), which can be explained through the DTS and FK models.   [21]. Similarly,  cr 2 becomes length-independent as L is larger than the characteristic width of solitons and their formation can be accommodated. The underlying mechanism of elastic effects is thus clarified as that solitons are nucleated as shear localized at the ends reaches the critical value,  cr [38]. The solitons and dislocations start to grow as the local elastic strain energy is large enough to overcome the energy barrier to nucleate dislocations (    10 1 8 10 J/m for  0


, corresponding to tension solitons) [40]. For very thick mesa and substrate, the penalty of elastic deformation is high, and overall sliding in the rigid limit could be achieved as the integrity of graphite flakes are well preserved. However, one should be noted that the deformation could be localized in the graphene layers adjacent to the interface, and the definition of  D t fails in the analysis of friction forces.
The analysis above can be extended to non-Bernal contacts between graphite mesa and substrate. For the "incommensurate" state with a misfit angle of | https://mc03.manuscriptcentral.com/friction     21 8 , the interlayer cohesive energy is Å 2 6 meV , which is one order lower than the value of Å 2 44 meV in Bernal stacking (Fig. S1(b) in the ESM). The interlayer shear modulus is measured to be   I 0 1 GPa G from the simulations, and  cr1 is estimated to be 0.1% and decreases with N, which is lower than the value of 0.6% at commensurate interfaces. As a result, elongation in graphite flakes is negligible, and the distribution of 2D atomic stress  ( ) x x is relatively uniform, with a small amplitude of 0.134 N/m (Figs. 4(c) and 4(f)). Moreover, corrugation in the interlayer potential energy surface is below . The width of strained soliton is thus predicted to be over 1,000 nm, much larger than the value of commensurate interfaces, and  cr 2 is less than 0.14%. The whole interface features a high level of lattice mismatch, and mechanical energy dissipation through lattice excitation is negligible [37]. The state of structural superlubricity can thus be preserved to a much extended length scale [2].
To validate the arguments, a set of MD simulations are conducted for interlayer sliding over the mesasubstrate contact with     21 8 (Fig. 4), which shows no reconstruction of the moiré pattern. For order lower than that of Bernal interfaces, arising mainly from the ends of contact [25,41]. It should also be noted that this value is of the same order as, although smaller than that measured from a rigid model (    3 6 5 10 GPa). No nucleation of local Bernal stacking is identified in the simulations. The negligible elongation of graphite flakes (less than 0.15%) is attributed to thermal activation processes during the loading process, which may diminish at low loading rates or in the quasi-static condition [42]. The out-ofplane displacement is also dominated by the end effects, with an amplitude of Å  0 3 , much more significant than that contributed by the moiré pattern ( Å  0 01 ). For thicker graphite (  5 N ), the friction stress is almost the same. However, the standard deviation, which measures the shear resistance, increases to    3 2 6 10 GPa, which may be attributed to the enhanced stiffness of mesa that restricts lattice relaxation. The lattice deformation is further suppressed to be under 0.05%, and the out-of-plane displacement is reduced to about 0.001 Å . Therefore, the rigid assumption applies for graphite contact with incommensurate interfaces. Increasing the thickness of graphene flakes elevates the in-plane stiffness, and further promotes the robustness of structural superlubricity [43].

Discussion
Elasticity at non-Bernal and incommensurate interfaces may lead to structural relaxation of the moiré pattern, resulting in the formation of local Bernal regions. This effect was reported for twisted bilayer graphene at "magic" stacking angles (    1 1 ), demonstrating intrinsic unconventional superconductivity, as well as renormalized vibrational and electronic structures [44,45]. Recent studies identified local atomic reconstruction up to a crossover angle of rigid twist at 5°, which results in abnormal electrical conductivity [46]. Our simulation results for a bilayer with     0 7 show that the moiré pattern selforganizes into a reconstructed triangular network by sacrificing the elastic deformation in graphene layers. The local atomic strain is 0.8%, and the out-of-plane displacement is Å    0 4 z and localized at ABB region. The reconstructed features remain intact during sliding as lattice distortion in the mesa is partially constrained by the loading condition (Fig. 5). The friction stress is fluctuating in a range of   0 2 GPa, much larger than the value measured at     21 8 for the appearance of Bernal-stacking regions, and approaching the value for Bernal interfaces (0.18 GPa for sliding along the AB-SP-BA path in the rigid model). Relaxing the constraint by applying the displacement to a few of the carbon atoms in the mesa leads to changes in the moiré pattern, reducing the friction stress. The structural superlubricity is thus perturbed at small twist angles as lattice reconstruction is activated.
Our atomistic simulation results show that elastic relaxation reduces the friction stress for Bernal contacts, but has negligible effects on the non-Bernal and incommensurate ones (Fig. S8 in the ESM). In continuum models, friction between elastic plates or at tip-substrate contacts is usually modeled by constitutive rules (e.g., the Coulomb law) implemented through interaction between surface nodes with chosen damping coefficients [47][48][49], or the cohesive zone model (CZM) with parameters derived from interatomic potential functions [50]. However, evolution of the lattice registry has not been incorporated in these models, and further work is needed to capture the multiscale features we report here. Crystal plasticity models with parameters controlled by the dislocation dynamics or evolutionary dynamics of strained solitons could be one of the choices [51].
In addition to the transverse force measured against relative sliding between contact pairs in simulations, friction can alternatively be quantified through the process of mechanical energy dissipation. Rigid models introduce mechanical energy dissipation of the system through phenomenological damping terms. Moreover, in reduced models such as PT and FK, the coarsegrained representation of lattices results in incomplete sampling of the phase space and thus cannot capture the essential process of dissipation. It should be noted, although lattice distortion reduces friction stress at interfaces with initial Bernal stacking, the lattice slips and phononic excitation play an important role in the kinetic friction of graphitic superlubricity systems [10]. At high sliding speeds, catastrophic damping effects were reported in the resonant regime [37,52,53]. With the size of contact much beyond the width of strained solitons, the nucleation of domains with local Bernal stacking could increase mechanical energy  | https://mc03.manuscriptcentral.com/friction dissipation through this regime. Quantitative discussion on this process, however, awaits further work.
It should also be remarked here that the simulations and theoretical models in our work and the literature are mostly restricted by their one-dimensional (1D) setup. With the three-dimensional (3D) nature of graphitic contacts considered, additional effects on the coupling between lattice slip and deformation have to be included in the discussion. The first one is the extended 2D nature of lattice slips and strained solitons [54]. The second effect is the rotation of mesa over the substrate [55]. The side effects for the finite-size mesa could also be prominent [16]. These complexities require more intensive theoretical and computational studies, as well as further exploration of model systems of structural superlubricity in experiments.

Conclusions
The limit of structural superlubricity can be achieved by breaking the interlayer lattice registry. The effects of elastic deformation and induced changes in the lattice registry are discussed here by exploring the nucleation and development of out-of-plane dislocations at the interface as well as strained solitons in the graphene layers. The size of contact and thickness are discussed as two controlling parameters that modulate the elastic deformation in graphene layers. By studying the sliding behaviors between graphite flakes using molecular dynamics simulations, we find that with an initial state of Bernal stacking, dislocations develop at the interface and even between the layers next to the interface, leading to significant elastic deformation. At non-Bernal interfaces with much larger superlattice periodicities or incommensurate contacts, the lattice distortion during sliding is negligible. To accommodate the formation of strained solitons or out-of-plane dislocations, a much larger contact than a micron is required. It is also interesting to note that, even with the formation of local commensurate regions, the friction stress with initial Bernal stacking is only one order higher than the non-Bernal contact, and reduced by one order from the value under the rigid assumption. For interfaces with very small twists, lattice selforganizes at the interface, and the reconstructed moiré pattern elevates the friction stress. These results thus approve the rigid assumption used in previous discussions on graphitic superlubricity at microscale [2,3]. The robustness of structural superlubricity can be improved for contact pairs with a small size of contact or very large stiffness of thick flakes. The work also adds insights into theoretical modeling of structural superlubricity, and formalizes the design recipes of practical applications.