Linking Friction Scales from Nano to Macro via Avalanches

Steady-state fluctuations in the friction force of molybdenum disulfide (MoS2), a prototypical lamellar solid, were analyzed experimentally for newton-scale forces and computationally via molecular dynamics simulations for nanonewton-scale forces. A mean field model links the statics and the dynamics of the friction behavior across these eight orders of magnitude in friction force and six orders of magnitude in friction force fluctuations (i.e., avalanches). Both the statistics and dynamics of the avalanches match model predictions, indicating that friction can be characterized as a series of avalanches with properties that are predictable over a wide range of scales.

Understanding the shear strength of lamellar solids like MoS 2 and graphene, materials that are widely used for their exceptional lubricity (i.e., low shear strength), is important for many engineering applications where traditional fluid lubricants cannot be used. Commonly used in mechanisms that operate in inert environments, such as in hermetic enclosures or exposed to the vacuum of space, MoS 2 films in the form of engineered coatings can undergo a transition to low friction coefficients (µ) during use, a phenomenon commonly referred to as "run-in." This transition in friction behavior is accompanied by structural changes of the coating surface [1][2][3] and the formation of a transfer film [4][5][6] on the mating counterface, generating a self-mated contact. The first few cycles of sliding cause shear-induced reorientation of the constituent atoms, leading to the formation of large, basally oriented MoS 2 sheets at the shearing interface parallel to the direction of sliding. These self-mated sulfurterminated lamellae interact through weak van der Waals forces that easily shear and provide low steady-state friction coefficients (µ < 0.05). In cases where the MoS 2 films are sufficiently pure, devoid of adsorbate contaminants, and basally oriented, superlubricity can be achieved, which is associated with friction coefficients µ < 0.01 [7]. Additionally, in this low friction regime, wear rates are also low, corresponding to the removal of as little as one molecular layer per sliding pass [8][9][10]. When in the steady-state friction regime, shearing planes of highly ordered MoS 2 can display persistent fluctuations in friction force (as shown in Fig. 1a and b) [8,11] that have been attributed to stick-slip behavior [12,13]. In this regime, atoms in adjacent lamella periodically become stuck in mutual local potential wells.
Here, we try to understand these low friction coefficients and wear rates from the statistics of fluctuations in the friction force that arise from underlying slip avalanches. An avalanche describes the cascading nature of failure events, experimentally observable as sharp drops in friction force. To understand the statistics of these avalanches, we model friction of MoS 2 lamellae via the slipping of elastically coupled weak spots between material layers [12,14,15]. These weak spots arise from the natural variations at the interface between the two sliding layers. Weak spots are regions that are most prone to slipping. Without thresholds to sliding, sliding would proceed smoothly in response to an applied external force. Heterogeneities at the interface, for example via local asperities, introduce a local threshold force to sliding, implying that the interface will remain stuck (or 'pinned') for small external forces. In response to a small force increase the interface slips in avalanches until it is pinned again by new asperities. Only for forces larger than a critical force, i.e., large enough to overcome all asperities, will the interface keep slipping. That critical force constitutes a critical point [16] that separates the pinned regime from the sliding regime. If we drive the system with a soft loading spring the system self-organizes to operate near the critical force. In that case it responds with avalanches that are widely distributed in size with statistical distributions that are characterized by scaling exponents and scaling functions. These only depend on a few basic properties such as symmetries, dimensions, interaction ranges, etc., but not on the microscopic details of the system. Below we describe a simple microscopic model that captures the basic properties and makes predictions for these universal exponents and scaling functions that can then be compared to experiments.
The statistics of these friction force avalanches resemble those in a chaotic regime of friction described in [10], as well as sudden stress drops arising from internal friction during deformation of a wide range of solids across a broad range of length scales [12,[17][18][19][20][21][22]. The model suggests that this behavior is similar because the underlying mechanism in all cases is slip avalanche weak spots. Predictions from this model have been confirmed by recent experiments on internal friction during slow compression of bulk metallic glass specimens; in these materials, the weak spots are shear transformation zones (STZs) [16] that predominantly slip in shear bands. Here, we show the statistics and dynamics of friction force avalanches measured in macroscale tribological experiments on sputter-deposited MoS 2 thin films. The results are compared with data from nanoscale molecular dynamics (MD) simulations of shear between MoS 2 lamella and predictions from a simple analytic model for the slip statistics during friction. The application of renormalization group methods to this simple analytic model enables the determination of aspects of the slip statistics and dynamics that are Friction forces from MD simulations show similar fluctuations. c The force-drop rate profile of the single experimental avalanche and the maximum force-drop rate Γ max . We use only experimental data from the steady-state regime, as shown in (d). Note in (b), the flat sections are not actually flat, but slightly curved downward, and thus the region between the red lines is one single avalanche independent of microscopic details and therefore apply across a wide range of scales. We find good agreement between the experiments, simulations, and theoretical mean field theory (MFT) model.
In plastic deformation experiments, avalanches manifest as sudden drops in stress [12,[18][19][20][21][22] with statistics and dynamics that are described well by a simple MFT model of coupled weak spots undergoing slip [12]. For friction, avalanches are observed as sudden drops in friction force with weak spots corresponding to local asperities at the interface. We define u(r, t) (sometimes referred to as the displacement field) as the interfacial offset since time t = 0 of a weak spot at location r. With an applied shearing force F, elastic forces from all other asperities int (r, t), and asperity pinning due to pinning forces f R u, r, u r, t ′ < t , the displacement field u(r,t) evolves in time according to [12] where η is an effective viscosity-like constant that sets the time scale of the slips and depends on the microscopic details of the system. Renormalization group calculations show that elastic interactions are sufficiently long range so that mean field theory, which approximates int (r, t) to be a constant, gives the correct scaling behavior in the physical dimension for slips along a shear band or frictional surface [23]. Thus, we take the mean field limit by assuming each site couples to the average force of all other sites , where J characterizes how strongly the slipping weak spots interact, i.e., it really is an "elastic coupling coefficient." To make contact with the physical system of interest, the change in local friction of a weak spot slipping at position r′ is defined as where K represents the effective bulk loading stiffness. The friction force at any asperity f (r � , t) can surpass a failure threshold value f fail and slip to reduce its friction force to an arrest force f arr (r � , t) . The elastic interactions mediate the transfer of this released energy to other asperities and may trigger more slips in what we refer to as an avalanche. To mimic periodic noise during an experiment (e.g., from noise in the electrical circuits to measure forces or ringing of the force sensors), we overlay friction force traces produced by our model with a signal resembling the periodic noise. Friction force fluctuations of MoS 2 were measured in ultra-high vacuum (UHV; 1 × 10 −9 torr) with a linearreciprocating micro-tribometer at a 10 kHz acquisition rate [24] (Fig. 1a and c). Force measurements were achieved via capacitance probes that measure the deflection of a custom dual-axis cantilever. On the end of the cantilever, a spherical probe (the counterface) is held fixed while the substrate is displaced parallel to the contact plane using a motorized stage. To establish a reference dataset for system noise, measurements were performed while the probe was out of contact and the stage was moving, and also while in contact with the stage stationary. More information regarding measurement uncertainty for force measurements can be found in the previous manuscripts [25]. Experiments consisted of bidirectional sliding cycles at a normal force of 0.2 N and sliding speed of 1 mm/s. The 3.2-mm-diameter spherical probes, made of 440C stainless steel, were brought into sliding contact against 440C stainless steel polished coupons with an approximately 100-nm-thick pure MoS 2 film deposited using an N 2 spray impingement process [11]. Each experiment consisted of 500 bidirectional cycles with a 2-mm-long stroke length. MoS 2 coatings under shear typically show a run-in period where friction coefficients and wear rates evolve from high (µ > 0.1) to low (µ < 0.05) as the surface structure transitions from amorphous (or near-amorphous, depending on the deposition method) to ordered lamellae [26]. To capture the physics of the steady-state regime, only the final 400 cycles of each experiment were used to extract magnitudes S, lengths T, and force-drop rate Γ of the friction force avalanches (Fig. 1d). The method to extract the avalanche sizes and durations from the data is as follows: zero crossings of the time derivative of the friction signal ( Γ(t) ) from positive to negative indicate when avalanches (i.e., force drops) start (t = t start ). At the following zero crossing from negative to positive the avalanche ends (t = t end ), see Fig. 1a, b. The avalanche duration T is then defined as T ≡ t end − t start , and the avalanche size S is equal to the amount by which the friction force dropped during the duration, S ≡ F t start − F(t end ). Friction force microscopy (FFM) measurements were performed on an Asylum MFP-3D atomic force microscope (AFM). MoS 2 -coated silicon substrates were mounted inside a humidity sensing cell on the sample stage. The humidity cell was attached to high-purity nitrogen (N 2 ), such that the cell was maintained at a positive gas pressure and the relative humidity was close to that of the dry N 2 source. Rectangular quartz-like cantilevers with SiO 2 colloidal spheres were used to facilitate a well-defined contact geometry (Nanosensors CP-qp-SCONT-SiO probes with a sphere radius of R = 1 μm). The normal force calibration C z and spring constant k z for a single probe were measured with the thermal fluctuation method [27]; multiple C z and k z values resulted in means and standard deviations of C z = 1.48 ± 0.02 nN/V and k z = 14.6 ± 0.2 pN/nm. The lateral sensitivity S of the probe was assessed with the extended wedge method [28]; multiple S values resulted in S = 12.8 ± 1.9. Following calibration, the probe was brought into contact with the MoS 2 until a normal photodiode voltage V z = 5.26 V was reached, translating to a normal force F n = 7.8 nN. Repetitive, single-line scans were performed over a 1 μm distance at a frequency of 0.10 Hz. The lateral photodiode voltage during the trace V xt and retrace V xr were measured for each cycle. For each V xt and V xr pair, the friction loop width W = (V xt − V xr )/2 and offset Δ = (V xt + V xr )/2 were calculated to determine friction force F f = C z (S)(W) and coefficient of friction μ = (S)(W)/V z [28]. The F f data from the steady-state μ region (final 400 cycles) were analyzed for avalanches.
MD simulations were performed using a reactive force field [29] at a fixed shear velocity of 1 m/s. Two sets of simulations were performed: (1) one with MoS 2 single crystals (15 square lamellae of MoS 2 ≈ 9 nm on a side) to study commensurate (high friction) contact and (2) another with a single sheet rotated by 30° to study incommensurate (low friction) contact; details of the procedure are in [30]. Friction forces from the simulations showed persistent fluctuations like those in the experiments (Fig. 1b). This enabled a determination of sizes, durations, and force-drop rates for the avalanches in the steady-state sliding regimes (after ≈ 100 ps and ≈ 35 ps of sliding for the commensurate and incommensurate cases, respectively).
The distribution P(S) of force-drop sizes S (Fig. 1a, b) scales as P(S) ~ S −τ , with universal (i.e., detail-independent) exponent τ. The prediction for τ depends on whether there is slip weakening (τ = 2) or not (τ = 3/2). This scaling holds for S within a power law scaling interval S ∈ [S min , S max ]. The interval is bounded by an upper cutoff S max that depends on experimental tuning parameters such as static and dynamic friction, machine stiffness, sample size, and driving rate. The region [S min , S max ] is often referred to as the scaling regime, and the exponent τ is a critical exponent.  15 s, respectively, and the maximum force-drop rates by G 0 = 0.54, 4.64 × 10 -8 N/s, 1.03 × 10 -7 N/s, and 1.5 × 10 -4 N/s, respectively. We extracted several universal scaling properties from this regime and compared them to the critical exponents and scaling functions predicted by the MFT model (Fig. 2).
To avoid ambiguities from varying histogram bin widths [31], we focus on power laws in the complement of the cumulative distribution function (CCDF), or one minus the integral of the probability distribution function (PDF) P(S) . The CCDF of avalanche size C(S) is the probability to have an avalanche of size S or larger, and scales as C(S) ∼ S −( −1) (Fig. 2a) with τ = 3/2 in the absence of slip weakening. Similarly, the CCDF of avalanche durations C(T) scales as C(T) ∼ T −( −1) , where T is the bin size (Fig. 2b) with α = 2. Avalanches are predicted to proceed in time such that the distribution of sizes scales as ⟨S⟩(T) ∼ T 1∕ z (Fig. 2c), where angle brackets denote the average.
The size-averaged maximum force-drop rate scales as ⟨Γ max ⟩(S) ∼ S z (Fig. 2d). Historically, 1 z is the fractal dimension of the avalanches in the scaling regime, with σ describing the stress dependence of the largest avalanche in the scaling regime, z describing the scaling of the duration of avalanches with their diameter, and describing the scaling of the correlation length near criticality. Within error bars, all critical exponents obtained in experiments and simulations were consistent with the MFT simulations with emulated machine noise (see Supplemental Information, section c). To strengthen agreement, we compared to stricter predictions. Two of the most stringent validations are the time profiles of the average force-drop rate as a function of size and duration, or Γ-profiles [12,16,32]. The Γ-profiles of N = 21 avalanches within a bin size centered at S i were averaged to obtain the size-averaged Γ-profile is the force-drop rate of the nth avalanche in the ith bin, where N is the number of avalanches in the bin S i . Time was offset such that all avalanches start at t = 0 (t j = t − t 0 , where t 0 is the avalanche start time). The MFT model predicts that the averaged force-drop rates, drawn from different size bins in the scaling regime of C(S), have the same form as the universal scaling function [12,14] Figure 3a, c, and e (left column) show that plotting collapsed the profiles. The experimental noise caused deviation from the analytical prediction at small values of t/S 0.5 , however, when the experimentally determined periodic noise signal was included in the MFT model, the deviations were nearly the same (Fig. 3c and Supplemental Information). Matching the functional form of scaling functions was a more stringent test of the theory than merely matching power law exponents. As shown below, we find good agreement between experiments and theory for two different scaling functions and four exponents.
In addition, C(T) followed a power law in the duration scaling regime (Fig. 2b). Averaging the force-drop rate time profiles of avalanches of about the same duration T i gave the duration averaged Γ-profiles ⟨Γ T i � t j � ⟩. MFT predicts that avalanches in this regime obey the scaling law Figure 3b, d, and f (right column) show that this prediction also matches experiments and MFT with added noise. The noise in the experiments again obscured the shapes, here for both small and large t/S 0.5 , and the MFT with noise showed similar deviations (Fig. 3d). In summary, we see strong agreement between the experiments and the model predictions for (1) the avalanche statistics, as seen in the scaling exponents of the statistical distributions, and (2) the avalanche dynamics, described by the average Γ-profiles that match the predicted critical scaling functions. This agreement between theory and experiments indicates that the fluctuations in the friction force indeed reflect the underlying critical point predicted by our simple MFT model. That is an important result because the model predicts that the statistics and dynamics of these fluctuations should be universal and independent of the microscopic details and the scale of the system. In other words, the MFT predicts that friction experiments on different materials should yield similar results, and that the results from the lab scales should be transferrable to systems on much larger or much smaller scales, as long as the symmetries, dimensions, and interaction ranges are the same.
System oscillations paired with periodic noise limited our analysis to regimes away from spikes in the power spectrum (see Supplemental Information). Although the scaling regime in the experiments spans less than two orders of magnitude, the agreement of our data with more than six measures including scaling collapses with universal scaling functions suggests that the data strongly support model predictions. It is the robust scale invariance predicted by the model that explains why both UHV experiments at mN scales and MD simulations at nN scales follow the same MFT scaling behavior. As predicted by MFT, both experiments and simulations showed a distinct scale invariance, reflecting an underlying critical point. The agreement in The MFT model is not necessarily expected to describe experimental results for small event sizes and durations, as it is designed to describe only long wavelength and low frequency regimes where universal behavior is found. At small scales, microscopic details come into play and deviations between MFT and experiments are unsurprising. Such deviations can be seen in CCDFs (Fig. 2a, b) for small sizes and durations. Interestingly, in the small size and duration regime, experimental results more closely matched MD simulations with incommensurate contacts. While the registry of experimental systems is unknown, the two simulation cases (i.e., commensurate and incommensurate) represent bounding cases and suggest that macroscale steady-state low friction consists of primarily incommensurate shear (see Supplemental Information).
Moreover, there was no evidence for critical avalanches in the AFM FFM data (Fig. 2). However, this was again expected, considering the single-asperity nature of FFM [6,8,9]: Avalanches in FFM involve slipping of the entire contact area, i.e., the entire tip slides on the surface until it is caught by a strong enough asperity. The model presented here predicts that power law scaling is only seen when the slipping size is much smaller than the size of the contact area. If that is not the case then the model predicts that finite size effects prevent the emergence of power law scaling, because the tip size sets a length scale that destroys scale invariant behavior. We estimated a ratio of slipping areato-contact area, for the largest avalanche, and found it to be of order 1000 to 1 for the AFM experiments. Thus, these dramatic finite size effects explain the absence of scaling of avalanches in AFM data. On the other hand, for the MD simulations, we find this ratio to be 1.0 and for the UHV experiments it is 2.5. If the largest slip is comparable to the contact area, then all smaller slips in the experiment are smaller than the contact area, and power law scaling can emerge as predicted. We note that the MD simulations do not have multiple, simultaneous slip locations, but instead the whole interface slides at once [33]. Importantly, there are layers beyond the interface that also slip, so in this case our model can be understood as describing the slipping of elastically coupled layers.
Importantly, we point out that the AFM data can still be described by our model, albeit in a parameter regime where the slips are large compared to the contact area (or system size). Similar effects can be achieved in the model, for example for large slip weakening in the system, and in a somewhat different way also for effectively soft loading springs. Either way this type of behavior is seen in model parameter regimes where the system is not close to a critical point and thus the avalanches are not power law distributed. Thus, these model predictions match the observations from the AFM experiments also.
Finally, new modeling work (arxiv:2204.02795) further supports the idea that friction force fluctuations should match our model with weakening [23,34], based on comparison to simulations of a rate and state model [35]. We note importantly here that, though the slope of the size and duration CCDF in the mentioned work are steeper than the MFT prediction without weakening (instead being a closer match to the prediction of the model with weakening), our MFT + noise simulations suggest this is simply a consequence of the periodic machine noise. It will be interesting to see if the exponents resembling those of the weakened model will persist in experiments without machine noise.
In summary, the results demonstrated that friction behavior of lamellar solids, and possibly materials in general, can be accurately described by our MFT model, and that this framework can be used to understand the similarities and differences between macroscale and nanoscale friction measurements. The framework can also be used to transfer friction results between scales, predict future friction behavior, and diagnose and quantify potential changes in friction over time. As friction is linked to plastic deformation, the broad agreement with predictions from the MFT model suggests the possibility of linking friction force statistics with internal friction and wear rates as a future application.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s11249-022-01619-x. otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.