Stochastic Properties of Static Friction

The onset of frictional motion is mediated by rupture-like slip fronts, which nucleate locally and propagate eventually along the entire interface causing global sliding. The static friction coefficient is a macroscopic measure of the applied force at this particular instant when the frictional interface loses stability. However, experimental studies are known to present important scatter in the measurement of static friction; the origin of which remains unexplained. Here, we study the nucleation of local slip at interfaces with slip-weakening friction of random strength and analyze the resulting variability in the measured global strength. Using numerical simulations that solve the elastodynamic equations, we observe that multiple slip patches nucleate simultaneously, many of which are stable and grow only slowly, but one reaches a critical length and starts propagating dynamically. We show that a theoretical criterion based on a static equilibrium solution predicts quantitatively well the onset of frictional sliding. We develop a Monte-Carlo model by adapting the theoretical criterion and pre-computing modal convolution terms, which enables us to run efficiently a large number of samples and to study variability in global strength distribution caused by the stochastic properties of local frictional strength. The results demonstrate that an increasing spatial correlation length on the interface, representing geometric imperfections and roughness, causes lower global static friction. Conversely, smaller correlation length increases the macroscopic strength while its variability decreases. We further show that randomness in local friction properties is insufficient for the existence of systematic precursory slip events. Random or systematic non-uniformity in the driving force, such as potential energy or stress drop, is required for arrested slip fronts. Our model and observations...


Introduction
Static friction is the maximal shear load that can be applied to an interface between two solids before they start to slide over each other. The famous Coulomb friction law (Amontons, 1699;Coulomb, 1785;Popova and Popov, 2015) states that static friction is proportional to the normal load with the friction coefficient being the proportionality factor. The friction coefficient is generally reported as function of the contacting material pair, which is often misinterpreted as the friction coefficient being a material (pair) property. While proportionality of friction to normal load is mostly valid, the friction coefficient is geometry-dependent and thus varies for different experimental setups with the same material pair (Ben-David and Fineberg, 2011). The underlying cause for this observation is the mechanism governing the onset of frictional sliding, which has been shown to be a fracture-like phenomenon (Svetlizky and Fineberg, 2014;Svetlizky et al., 2020;Rubino et al., 2017). The geometry and deformability of the solids lead to a non-uniform stress state along the interface. As a consequence, local frictional strength is reached at a critical point and slip nucleation starts, from where it extends in the space-time domain -just like a crack -until the entire interface transitioned and global sliding occurs.
Variations in the static friction force, however, do not only occur because of changes in the loading configuration. Experiments have shown that the measured friction force varies also from one experiment to another when the exact same setup and exact same specimens are used. For instance, Rabinowicz (1992) showed that the static friction coefficient of a gold-gold interface, measured by a tilting plane friction apparatus, varies from 0.32 to 0.80 for normal load 75 g. Similar but to a lesser extent, Ben-David and Fineberg (2011) also observed variations in the static friction coefficient of glassy polymers when the loading configuration was fixed.
While these variations are not often reported, they are an important factor in the absence of a complete and consistent theory for friction (Spencer and Tysoe, 2015). If (seemingly) equivalent experiments lead to a large range of observations without consistent trends, it is challenging to isolate the relevant from the irrelevant contributions and, therefore, nearly impossible to create a fundamental understanding of the underlying process. Even though the presence of these large variations has important implications for the study of friction, current knowledge about the origin and properties of these observed variations in macroscopic friction remains limited.
One possible origin is randomness in local friction properties. Interfaces have been shown to consist of an ensemble of discrete micro-contacts (Bowden and Tabor, 1950;Dieterich and Kilgore, 1994;Sahli et al., 2018), which are created by surface roughness (Thomas, 1999;Hinkle et al., 2020) when two solids are brought into contact. This naturally leads to a system with random character, where micro-contacts of random size are distributed randomly along the interface (Greenwood et al., 1966;Persson, 2001;Hyun and Robbins, 2007;Yastrebov et al., 2015). Since frictional strength is directly related to the cumulative contact area of theses micro-contacts (Bowden and Tabor, 1950;Greenwood et al., 1966), and the micro-contacts are the result of random surface roughness, the local frictional strength is likely also random.
Surprisingly, the effect of interfacial randomness on friction remains largely unexplored. Most of previous work is focused on how (random) surface roughness is related to various friction phenomenology including rate-dependence (Li et al., 2013;Lyashenko et al., 2013), local pressure excursions within lubricated contact (Savio et al., 2016), chemical aging (Li et al., 2018), or the existence of static friction (Sokoloff, 2001). How interfacial randomness causes variation in these observations has, however, not been studied so far. Only most recently, Amon et al. (2017) and Geus et al. (2019) have considered variability in friction. Amon et al. (2017) showed that systems with a nonuniform initial stress state with long range coupling are characterized by two regimes: at low loading, small patches of the system undergo sliding in an uncorrelated fashion; at higher loading, instabilities occur at regular intervals over patches of increasing size -just like confined stick-slip events (Kammer et al., 2015;Bayart et al., 2016) -and eventually span the whole system. Geus et al. (2019) simulated interface asperities as an elasto-plastic continuum with randomness in its potential energy and show that the stress drop during a stick-slip cycle is a stochastic property which vanishes with increasing number of asperities. These results demonstrate well the stochastic character of macroscopic friction due to random interface properties. However, the effects of interfacial randomness on the variability of macroscopic static friction, e.g., the friction coefficient, has not been studied yet.
Here, we address this gap of knowledge and aim at a better understanding of the stochastic properties of static friction. We present a combined numerical and theoretical study that links randomness of local friction properties with observed variability in macroscopic strength. Using dynamic simulations, we will show that the macroscopic friction threshold is attained when a local slipping area, of which many can co-exist, reaches a critical length and nucleates the onset of friction. This nucleation patch becomes unstable and propagates across the entire interface causing global sliding. We will then show that a quasi-static equilibrium theory, which takes an integral form, predicts quantitatively well the critical stress level that causes nucleation of global sliding. Based on this theoretical model, we will develop fast and accurate Monte Carlo simulations using a Fourier representation of the integral equations, and demonstrate the extent of variability in macroscopic static friction based on random interface strength with various correlation lengths. Finally, we will show that a decreasing interfacial correlation length leads to higher macroscopic strength with decreased variability. This paper is structured as follows. First, we provide a problem statement in Sec. 2 including a description of the physical system, the stochastic properties, and our approach to generate random strength fields. In Sec. 3, we present the numerical method used to simulate the onset of frictional sliding and compare simulation results of critical stress leading to global sliding with predictions based on a theoretical model. This model is then used in an analytical Monte Carlo study, which is developed and presented in Sec. 4. The implications of our model assumptions as well as the model results are discussed in Sec. 5. Finally, we provide a conclusion in Sec. 6.  ) of strength τ f (δ, x) embedded within two semi-infinite elastic solids, which are periodic in x with period L and infinite in y. A uniform loading τ 0 (t) is applied. (b) τ 0 (t) is increased linearly with time t up to the onset of frictional motion when the stress drops from its critical value τ cr to a kinetic level τ kin . (c) The constitutive relation of the frictional interface is a linear slip-weakening law τ f (δ) with random peak strength τ p (x) and constant weakening rate W (see Eq. 1). (d) τ p (x) is a random field with spatial correlation C(ξ) (inset) and probability density function f (τ p ) (right).

Problem Statement
In this section, we first provide a description of the physical problem that we consider throughout this paper. We then describe the stochastic properties of the strength profile along the interface and, finally, explain how we generate these random fields.

Physical Problem
We study the macroscopic strength of a frictional interface. Our objective is to provide a fundamental understanding of the effect of local variations in frictional strength on the macroscopic response. For this reason, we focus on the simplest possible problem -without oversimplifying the constitutive relations of the bulk and the interface. Specifically, we consider a two-dimensional (2D) system consisting of two semi-infinite elastic solids, as shown in Fig. 1a. The domain is infinite in the y direction and periodic in x with period L. Both materials have the same elastic properties.
We apply a uniform shear load τ 0 (t) that increases quasi-statically with time (see Fig. 1b). Once τ 0 (t) reaches the macroscopic strength of the interface τ cr , the interface starts to slide and the frictional strength suddenly reduces to its kinetic level τ kin . This observed reduction in shear stress is typically associated with friction-weakening processes, which may depend on various properties, such as slip, slip rate, and interface state. The critical shear stress τ cr , if divided by the contact pressure, corresponds to the static macroscopic friction coefficient. Similarly, τ kin is proportional to the kinetic friction coefficient.
These macroscopic observations depend on the local interface properties, which are the peak strength τ p (x), and residual strength τ r (x). As we will show, the local properties are generally different from the macroscopic properties; particularly, in the case of non-uniform stress or strength. For simplicity, we describe the evolution of the local frictional strength as a linear slip-weakening law, which is shown in Fig. 1c and is given by where δ(x) is local slip, d c (x) is a characteristic length scale, and W(x) = (τ p (x) − τ r (x))/d c (x) is the weakening rate. H(.) is the Heaviside function.
We consider a heterogeneous system with local peak strength τ p (x) being a random field, as further described in Sec. 2.2. To reduce complexity of the problem, we assume uniform residual strength 1 ∂ x τ r = 0 and uniform weakening rate ∂ x W = 0. The variation in local peak strength is thought to represent possible heterogeneity in the material, but also the effect of surface roughness, which leads to a real contact area that consists of an ensemble of discrete contact points with varying properties. The implications of this approach will be discussed in depth in Sec. 5.

Stochastic Properties of Frictional Interface
The local peak strength τ p (x) is modeled as a stationary non-Gaussian random field with specified cumulative distribution function F(τ p ) and corresponding probability density f (τ p ), as shown in Fig. 1d. The random field is defined by the nonlinear mapping where z(x) is Gaussian with zero mean and unit variance and Φ its cumulative distribution, depicted in Fig. 2a-left. F and Φ are monotonic by definition, so their inverse exist, which can be used to prove that P τ p (x) ≤ τ = F(τ). Further, with τ p being the local peak strength of the interface, it needs to satisfy some physical requirements. First, the peak strength is always higher than the residual strength, i.e., τ min p ≥ τ r . Second, it maximum value is limited by the material properties. For this reason, we require that τ p ∈ (τ min p , τ max p ), which we achieve by setting F(τ p ) as a Beta cumulative distribution function (see Fig. 2a-right).
The spatial evolution of z(x) is specified by its power spectral density g(k), which corresponds to the Fourier transform of the correlation function C z (ξ), i.e., where k is the angular wave number. We assume that z(x) has a power spectral density where λ is the cutoff frequency, above which the spectral density decays as a power law ∼ k −8 (see Fig. 2b). The correlation length ξ 0 is a measure of memory of the random field; the longer ξ 0 the longer the memory. ξ 0 is inversely proportional to λ, and we define 2 it as ξ 0 = 2π/λ. Since the correlation function of z is positive, C z (ξ) > 0, it is not greatly affected by the nonlinear mapping F −1 • Φ and C z (ξ) ≈ C τ p (ξ) (Grigoriu, 1995, p.48) and features such as ξ 0 are preserved (see Fig. 1.d inset). The assumption of using this specific spectral density and probability distribution are discussed in Sec. 5.3.

Random Field Samples Generation
The samples of random field τ p (x) are generated as follows. First, the Gaussian random field z(x) is generated using a spectral representation where A j and B j are independent Gaussian random variables with zero mean and unit variance and modal angular wave-number is k j = 2π j/L. The fundamental wavelength of the field 2π/k 1 = L is chosen such that it corresponds to the domain size L, which implies that z(x) is periodic over L, and so is τ p (x). The modal variance σ 2 j ∝ g(k j ) corresponds to the discrete spectral density, which is normalized to assure that z has unit variance 1 We use ∂ i as short notation for partial derivative with respect to i. The random variable τ p is generated by applying a nonlinear mapping F −1 • Φ(x) (Eq. 2) onto the Gaussian random variable z. Two colored dots inked by a line represent a (z, τ p ) pair with equal cumulative density F(τ p ) = Φ(z), which is the criterion imposed by the nonlinear mapping. (b) Normalized spectral density function g(k)/g(0) (Eq. 4).k is the truncation frequency and λ = 2π/ξ 0 is the cutoff frequency.
Due to the discrete representation of z(x), we apply a truncation frequency that is considerably larger than the cutoff frequencyk ≡ k J = 2.5λ. This ensures that most of the spectral power is preserved: Further increase ink would include additional high frequency modes but with negligible amplitudes. Finally, once z(x) has been generated, we apply the nonlinear mapping F −1 • Φ (Eq. 2 and visualized in Fig. 2a) and obtain the random field τ p (x). Fig. 1d shows a sample of τ p (x) generated using the described procedure with corresponding correlation function and probability density.

Dynamic Simulations
In the following, we will first present the numerical method and model setup applied in our simulations of the onset of friction. We then provide a theoretical model to describe the simulations and present a comparison between the theoretical predictions with the numerical results.

Numerical Method
We model the physical problem, as described in Sec. 2.1, with the Spectral Boundary Integral Method (SBIM) (Geubelle and Rice, 1995;Breitenfeld and Geubelle, 1998). This method solves efficiently and precisely the elastodynamic equations of each half space. The spectral formulation applied in SBIM naturally provides periodicity along the interface. The half spaces are perfectly elastic and we apply a shear modulus of G = 1 GPa, Poisson's ratio of ν = 0.33 and density ρ = 1170 kg/m 3 , and impose a plane-stress assumption. While we will report our results in adimensional quantities, we note that these parameters correspond to the static properties of glassy polymers, which have been widely used for friction experiments (Svetlizky and Fineberg, 2014;Rubino et al., 2017).
The interface between the two half spaces is coupled by a friction law as given by Eq. 1. The friction law corresponds essentially to a cohesive law, as known from fracture-mechanics simulations, but applied to the tangential direction. It describes the evolution of local strength as a function of slip. We apply peak strength τ p (x) as a random field, following the description provided in Sec. 2.2, and constant τ r . τ p (x) follows a Beta distribution with α = 1.5 and β = 3. We impose a maximum value for relative peak strength of max(τ p (x) − τ r ) = 1.66 MPa and minimum value of min(τ p (x) − τ r ) = 0.66 MPa. Therefore, the random field has a mean value of τ p − τ r = 1 MPa and standard deviation of 0.2 MPa. We further apply a constant slip-weakening rate of W = 0.5 TPa/m, which is representative for glassy polymers (Svetlizky et al., 2020). Finally, a slowly increasing uniform stress τ 0 (t) is applied along the entire interface. We use a repetition length of L = 0.1 m, which is, as we will show, considerably larger than the characteristic nucleation length scale. The interface is discretized by 512 − 1024 nodes. We verified convergence with respect to discretization, loading rate, and time step.
The results of a representative simulation are shown in Fig. 3. The τ p (x) profile has many local minima (see Fig. 3c). Depending on their value, these minima cause localized slip, as evidenced by bright blue vertical stripes over most of the time period shown in Fig. 3b-right. These localized slip patches grow slowly with increasing loading, which is difficult to see for most patches in Fig. 3b-right. Growth is easiest observed for the slip patch at x/L ≈ 0.75. Incidentally, this patch grows enough to reach a critical size from which on the patch becomes unstable, marked by a black dot, and starts growing dynamically. This dynamic propagation, see orange-red area in Fig. 3a enlarged from Fig. 3b-right, does not stop and, therefore, causes sliding along the entire interface -hence global sliding. The effect on the macroscopic applied force is shown in Fig. 3b-left, where τ 0 (t) = L τ(x, t)dx. At the precise moment when the slip patch becomes unstable, marked by a black dot, τ 0 (t) starts decreasing rapidly. The maximum value, denoted τ cr , represents the macroscopic strength of the interface.
The simulation shows that macroscopic strength is not reached when the first point along the interface starts sliding but when the most critical slip patch becomes unstable, starts propagating dynamically, and "breaks" the entire interface. Therefore, the criterion determining macroscopic strength is non-local and depends on the stability of local slip patches. In the following section, we will present a theoretical description of this nucleation process and provide a criterion for the limit of macroscopic strength.

Theory for Nucleation of Local Sliding
During the nucleation process, a weak point along the interface starts sliding. Due to local stress transfer, the size of this slipping area grows continuously until it reaches a critical size and unstable interface sliding occurs (Campillo and Ionescu, 1997). In this section, we will adapt the criterion developed by Uenishi and Rice (2003), which is shortly summarized in Appendix A, to describe and predict the limits of stable slip-area growth. Uenishi and Rice (2003) considered a similar system with two main differences to the problem studied here. First, in their case, the interface strength is uniform and the applied load is non-uniform. Appendix A shows that both problems result in the same equation for the problem statement and thus lead to the same nucleation criterion. Second, Uenishi and Rice (2003) considered a system with an isolated non-uniformity in the applied load. In other words, the applied stress was mostly uniform but with one well-contained local increase. Therefore, the location of nucleation is known in advanced. In our system, where the non-uniform property is random, the location is unknown. We will address this difference here and discuss it further in Sec. 5. Uenishi and Rice (2003) showed that on interfaces governed by linear slip-weakening friction (Eq.1), there is a unique critical length for stable growth of the slipping area, which can be approximated by where G * = G/(1 − ν) for mode II plane-stress ruptures, assuming the stress within h n has not attained the residual value τ r anywhere. Eq. 8 shows that h n depends only on the shear modulus G * and the slip-weakening rate W. Most importantly, the critical length is independent of the shape of the non-uniformity in the system. Specifically to our case, it does not depend on the functional form of τ p (x). Since we have homogeneous elastic solids and a uniform slip-weakening rate W, the critical size h n is unique and uniform along the entire interface. The important question for our problem, however, is to determine the level of critical stress that causes a nucleation patch to reach h n and initiate global sliding. The solution for the stress level leading to nucleation, as derived by Uenishi and Rice (2003), is given by Eq. A.7, and can be rewritten in terms of τ p (x) and h n as where x is the center location of the nucleation patch and v 0 (s) ≈ (0.925 − 0.308s 2 ) √ 1 − s 2 is the first eigenfunction of the elastic problem. Note that the transformation applied to the argument of τ p (x) results in the integral being computed over the critical nucleation patch size h n . Eq. 9 shows that the nucleation stress, which leads to a nucleation patch of size h n , does clearly depend on the shape of τ p (x). Note that Eq. 9 assumes that stresses within the nucleation patch have not attained the residual strength yet i.e., the fracture process zone spans the entire crack. In this case, the assumption of small-scale yielding does not hold, and, thus, the Griffith criterion for crack propagation does not apply.
As stated earlier, the nucleation stress τ n (x) was derived for a contained non-uniformity, for which we know the location. Therefore, τ n (x) corresponds to the critical stress of the system. In our system, however, τ p (x) is random and multiple nucleation patches might slowly grow. Determining the critical stress τ cr of the system requires computing the nucleation stress τ n for each nucleation patch and identifying the critical one. To address this aspect, we propose to compute Eq. 9 as a weighted moving average over the entire interface, and define the critical stress to be its minimum (see Fig. 4a). Therefore, we define the critical stress τ cr as For simplicity, we refer to this definition also as τ cr = min(τ n (x)) and x cr = arg min(τ n (x)). While it is possible, but not very likely, to have multiple minima of τ n with the same amplitude, this does not affect the resulting τ cr . However, multiple x cr could coexist which would result in multiple slip patches becoming unstable simultaneously. By adopting Eq. 9 and defining Eq. 10, we essentially assume that there is no interaction between nucleation patches. We will verify the validity of this assumption in the following section.

Results
We compare the results from numerical simulations, as described in Sec. 3.1, with the theoretical prediction from Sec. 3.2 by analyzing simulations with random τ p (x) generated using the method described in Sec. 2.3. For each of the three different correlation lengths ξ 0 /h n = 0.25, 0.5, and 2.0 we run 20 simulations. The system size is fixed and chosen such that it is considerably larger than the nucleation length, i.e., h n /L = 0.034.
A representative example is shown in Fig. 3. The size of h n /L is indicated in Fig. 3c and appears to provide a reasonable prediction for the nucleation patch size as observed in Fig. 3a. Further comparison is given in Fig. 4a. First, we illustrate the theoretical prediction. The nucleation stress τ n (x) (solid blue line) is computed from τ p (x) (gray dashed line) using Eq. 9, and τ cr is, according to Eq. 10, the minimum of τ n (x) (marked by black dot). We find the location of nucleation to be x cr /L ≈ 0.75, which corresponds to our observation from the numerical simulation, as seen in Fig. 3.
A more precise and systematic comparison is provided in Fig. 4b&c. We compare the predicted critical stress τ pred cr with the measured value from dynamic simulations τ sim cr . We compute τ pred cr as described above with Eq. 10, and as illustrated in Fig. 4a. We further find τ sim cr = T ∂ t τ 0 , where ∂ t τ 0 is the applied loading rate and T is the time at which τ 0 (t) is maximal (see Fig. 3b-left). Comparison of τ pred cr with τ sim cr is shown in Fig. 4b for all 60 simulations. The results show that the prediction works generally well. For decreasing ξ 0 /h n the prediction becomes slightly less accurate with a tendency to over-predict the critical value. The results further show that the predicted and measured critical stress τ cr increases with decreasing ξ 0 /h n .
While the location of nucleation is not relevant for the apparent global strength of our system, we compare the predicted and simulated x cr for further evaluation of the developed theory. The comparison shown in Fig. 4c uses x pred cr , as given by Eq. 10 and shown for an example in Fig. 4a, and x sim cr as found by analyzing the simulation data as illustrated in Fig. 3a&b-right. The data shows that the prediction works well for most of the simulations. For 8 simulations, 6 of which have ξ 0 /h n = 0.25, the prediction does not work. However, as shown in Fig. 4b, τ cr , which is the quantity of interest here, is correctly predicted for all of these cases. The reason for this discrepancies are likely second-order effects, as we will discuss in Sec. 5.
Overall, the results show that τ cr is quantitatively well predicted by the theory presented in Sec. 3.2. This allows us to study systematically the effect of randomness in interface properties by applying the theoretical model in analytical Monte Carlo simulations.

Analytical Monte Carlo Study
In the following section, we introduce Monte Carlo simulations, which are based on the theoretical framework for nucleation of frictional ruptures in a random field of frictional strength τ p (x), as derived in Sec. 3.2. The effect of correlation length ξ 0 on the effective frictional strength τ cr (Eq. 10), and its probability distribution f (τ cr ), is studied, while keeping all other properties constant. A Monte Carlo study based on the full dynamic problem (Sec. 3.1) would be computationally daunting. However, the theoretical framework allows us to evaluate τ cr very efficiently and has been validated by 20 full dynamic simulations for each considered ξ 0 (see Fig. 4).

Monte Carlo Methodology
The effective frictional strength τ cr = min(τ n (x)) requires the computation of the nucleation strength τ n (x), which involves a convolution of the local peak strength τ p (x) with the eigenfunction v 0 , given in Eq. 9. Considerable computation time can be saved by using a spectral representation of the random field τ p (x): where˜signifies thatτ p (x) is an approximation of τ p (x) and the number of frequencies J is chosen such that the approximation error |τ p − τ p | is negligible.τ p (k j ) is the discrete Fourier transform of τ p (x) where τ p (x) is generated using the procedure described in Sec. 2.3. By substituting Eq. 11 into Eq. 9 the nucleation strength convolution becomes a dot product: where g j (x) = +1 −1 e −ik j hn 2 s+x v 0 (s)ds is the modal convolution term, which, being independent of the sample specific functional form of τ p (x), can be pre-computed. This formulation allows for efficient and precise evaluation of the effective frictional strength τ cr = min τ n (x) for a large number of samples N = 10, 000, such that the probability distribution f (τ cr ) and its evolution as function of the correlation length ξ 0 can be accurately studied.

Monte Carlo Results
Prior to presenting the numerical results we provide some intuition of the effect of correlation length on the nucleation strength τ n based on probabilistic arguments. By exploiting the stationarity of τ p and τ n it is possible to derive an analytical expression of the expectation of the nucleation strength E[τ n ] and its variance Var[τ n ] as function of the corresponding statistical properties of local strength, E[τ p ], Var[τ p ] and ξ 0 /h n (see Appendix B).
One interesting finding is that the expectation is not affected by ξ 0 /h n : E[τ n ] = E[τ p ] (see derivation in Eq. B.2). The expression for Var[τ n ], however, involves a double integral of the product of the correlation function C(.) and the eigenfunction v 0 (.), which can be evaluated numerically (see derivation in Eq. B.3). For perfect correlation, i.e., ξ 0 /h n = ∞, C(.) becomes a constant, thus Var[τ n ] = Var[τ p ]. Additionally, in the limit of ξ 0 h n , the double integral in Eq. B.3 scales with ξ 0 /h n , thus Var[τ n ] ∝ Var[τ p ]ξ 0 /h n (see derivation in Eq. B.5).
We consider a range of correlation lengths ξ 0 /h n = {0.25, 0.5, 1.0, 2.0}, while all other properties remain constant. Fig. 5a-left shows one sample of τ p (x) for each considered ξ 0 . For clarity of visualization, in Fig. 5a we use the same seed when generating the random fields. Hence, the fields have the same modal random amplitudes A j and B j , see Eq. 5, but have different modal spectral densities σ 2 j , corresponding to the different ξ 0 . For this reason, all shown samples have a similar spatial evolution and the effect of varying ξ 0 can be clearly observed. By definition, all τ p (x) samples are drawn from the same probability distribution f (τ p ) (see Fig. 5a-center). Decreasing ξ 0 , moves the probability density of its minimum f (min(τ p )) towards the lower bound τ min p = 0.66 τ p (see Fig. 5a-right), because with lower correlation lengths it is more likely to visit a broad range of τ p values. Fig. 5b-left shows the corresponding nucleation strength τ n (x) for each of the local frictional strength fields τ p (x) presented in Fig. 5a, computed using Eq. 13. As mentioned before, τ n is essentially a weighted moving average of τ p with window size h n (see Eq. 9). Thus, most of the high frequency content of τ p disappears and the effect of ξ 0 on τ n is more subtle. One interesting feature is in the minima and maxima of τ n : increasing ξ 0 causes lower minima and higher maxima, because the moving average is effectively computed over an approximately constant field τ n ≈ τ p . Inversely, decreasing ξ 0 causes the opposite effect and τ n ≈ τ p .
This effect is more clearly visible by considering the distribution f (τ n ) shown in Fig. 5b-right. Increasing ξ 0 effectively puts more weight onto the tails of f (τ n ) (see ξ 0 /h n = 2.0 in Fig. 5b-right), and in the limiting case of ξ 0 /h n → ∞ the distribution of τ n will be the same as the one of τ p (analogous to Eq. B.4). On the other hand, decreasing ξ 0 puts weight on its mean τ p , making f (τ n ) similar to a Gaussian (see ξ 0 /h n = 0.25 in Fig. 5b-right) with variance proportional to ξ 0 (see Eq. B.5). In the limit ξ 0 /h n → 0 the distribution of τ n becomes a Dirac-δ centered at τ p . The described dependence of f (τ n ) on ξ 0 confirms the previously stated statistical arguments (see Appendix B for derivation).
Because f (τ p ) is skewed towards the lower bound of τ p so is f (τ n ); the larger ξ 0 the larger the skewness. For τ cr this effect is amplified by the fact that τ cr = min(τ n (x)) as depicted in Fig. 5c, causing τ cr to decrease with increasing ξ 0 . As noted in Sec. 3.3, the location where the critical instability occurs x cr is uniformly distributed over the entire domain as shown in Fig. 5d and is independent on ξ 0 .
We further analyze the effect of ξ 0 on the probability distribution of τ n and τ cr by reporting the mean, median and 25% percentile of the probability density function (see Fig. 6). We observe that the nucleation strength tends towards the mean peak strength for vanishing correlation length, lim ξ 0 /h n →0 τ n (x) = τ p , because the moving average in computing τ n (x) is evaluated over a window h n that appears infinite compared to ξ 0 (see Fig. 6a). Consequently, the effective strength also tends towards the mean of the local strength: lim ξ 0 /h n →0 τ cr = τ p (see Fig. 6b). Conversely, if ξ 0 h n , the moving average is computed over a window h n which vanishes, and thus (9) becomes the identity: lim ξ 0 /h n →∞ τ n (x) = τ p (x). In this case, the effective strength will be more likely to be close to the actual lower bound of the distribution lim ξ 0 /h n →∞ τ cr = min(τ p ) (see Fig. 6b). The transition between these two limiting cases is described by the results of the analytical Monte Carlo study, which are validated by 20 dynamic simulations for each ξ 0 /h n by reporting the mean effective friction strength (see inset in Fig. 6b). For ξ 0 /h n ≥ 0.5 simulations and theory coincide. However, for ξ 0 /h n = 0.25 the theoretical model slightly overestimates the effective friction strength. In Sec. 5, we will argue that the lower effective friction at small ξ 0 /h n is caused by interactions between neighbouring nucleation patches, which destabilize each other and lead to unstable growth at lower overall loading compared to the theoretical prediction.

Implications of the Physical Problem
The analyzed physical problem is simplistic and contains only the absolute minimum of a realistic system with a frictional interface -while still maintaining a rigorous representation of the constitutive relation of the bulk and the interface. The objective is to provide a fundamental understanding of the macroscopic effects on static friction caused by randomness in the local frictional properties. While many options exist to complexify the proposed system, we leave them for future work and focus here on the basics. Nevertheless, in this section, we will discuss some of these simplifications as well as their implications.
Randomness along the interface may have various origins including heterogeneity in bulk material properties and local environmental conditions (e.g., humidity and impurities). Prominent causes for randomness are geometric imperfections, which include non-flat interfaces and surface roughness. The real contact area, which is an ensemble of discrete micro-contacts (Bowden and Tabor, 1950;Dieterich and Kilgore, 1994;Li and Kim, 2008;Sahli et al., 2018) and is much smaller than the apparent contact area, introduces naturally randomness to the interface. Surface roughness is often modeled as self-affine fractals (Pei et al., 2005), which directly affects the size distribution of micro-contacts and local contact pressure. The resulting frictional properties are expected to vary similarly. This would typically lead to small areas of the interface with high frictional strength and most areas with no resistance against sliding, i.e., τ p = 0, since only the micro-contacts may transmit stresses across the interface. Therefore, at this length scale, one would expect the random strength field to be bound by zero at most locations, similar to the approach taken by Barras et al. (2019). However, in many engineering systems, the nucleation length is orders of magnitude larger than the characteristic length scales of the micro-contacts: nucleation lengths of ∼ 10 − 100 mm (Ben-David and Fineberg, 2011;Latour et al., 2013) and surface roughness lengths of ∼ 1 µm (Svetlizky and Fineberg, 2014). For this reason, we consider a continuum description with a somewhat larger length scale. In our approach, the frictional strength profile is continuous and varies due to randomness in the micro-contacts population without considering individual contact points.
Surface roughness and other local properties directly affect how frictional strength changes depending on slip δ, slip rate ∂ t δ, and state (Rabinowicz, 1995;Pilvelait et al., 2020). This is often modeled in phenomenological rate-and-state friction models (Dieterich, 1979;Ruina, 1983;Rice and Ruina, 1983). As discussed by Garagash and Germanovich (2012) and demonstrated by Rubin and Ampuero (2005) and Ampuero and Rubin (2008), the nucleation length scale of rate-and-state friction models approaches asymptotically the critical length h n used in this work and given by Eq. 8 if the rate-and-state friction parameters are favoring strong weakening with slip rate. However, if rateweakening becomes negligible, the nucleation criterion tends towards the Griffith's length (Andrews, 1976), which applies to ruptures with small-scale yielding. In this case, the frictional weakening process is contained in a small zone at the rupture tip and most of the rupture surface is at the residual stress level, which is different to the nucleation patches by Uenishi and Rice (2003), where the entire rupture surface is still weakening when the critical length is reached.
Since many engineering materials present relatively important slip-rate weakening friction, e.g., dynamic weakening of ∼ 1 MPa for glassy polymers at normal pressure of ∼ 5 MPa (Svetlizky et al., 2020) and, similarly, ∼ 1 MPa weakening for granite at normal pressure of ∼ 6 MPa , we considered a model system with strong frictional weakening. However, we neglect the complexity of rate-and-state friction, as extensively demonstrated by Viesca (2017, 2019), and apply a linear slip-weakening friction law at the interface because it has the most important features of friction, i.e., a weakening mechanism, while being simple and well-understood. The advantage is that the weakening-rate W is predefined. It further has a well-defined fracture energy Γ, which is the energy dissipated by the weakening process, i.e., the triangular area (τ p − τ r )d c /2 in Fig. 1c: Since W is constant in our system Γ varies with (τ p (x) − τ r ) 2 . A correlation between Γ(x) and τ p (x) can be expected given that any point that is stronger, i.e., increased τ p , is likely to dissipate more energy as well, i.e., increased Γ. Further, the linear slip-weakening law is contact-pressure independent, which may appear counter-intuitive based on Coulomb's well-known friction laws (Amontons, 1699;Coulomb, 1785). However, the contact pressure is, due to symmetry in similar-material interfaces, constant over time and, therefore, any possible pressure dependence becomes irrelevant for the nucleation process itself. Nevertheless, local friction properties are expected to change for systems with different normal pressure. This effect has not been analyzed here since we did not vary the contact pressure, but could be taken into account by changing the values of τ p , τ r and d c .
Finally, we note that by assuming a periodic system, we neglect possible boundary effects. We expect that the boundary would locally reduce τ cr compared to the prediction based on Eq. 10, which assumes an infinite domain, because the free boundary would locally restrict stress redistribution and thus increase the stress at the edge of the nucleation patch. Therefore, the probability density of global frictional strength τ cr for a periodic system, as shown in Fig. 5c, has likely a slight tendency towards higher values compared to a finite system. However, we expect the spatial range of the boundary effect to scale with h n and, therefore, f (τ cr ) will tend towards the periodic solution for h n /L → 0. Verification would require a large number of numerical simulations, which is beyond the scope of this work.

Interpretation of Numerical Simulations
The simulations have shown that uniform τ 0 and random τ p (x) cause multiple nucleation patches to develop simultaneously. We can see in Fig. 3b-right that 20-30 patches (bright blue stripes) coexist by the time global strength is reached, i.e., t/T = 1. Most of these nucleation patches grow very slowly and their number increases with increasing τ 0 (t). Nucleation patches can also merge, which is what happens in this simulation to the critical patch. Furthermore, the simulation shows that unstable growth and thus global failure is not necessarily caused by the first nucleation patch to appear. For instance, the τ p (x) profile shown in Fig. 3c presents three local minima with approximately the same value, i.e., at x/L ≈ 0.4, 0.6 and 0.75. Therefore, the first three nucleation patches appear quasi-simultaneously. Whether one of these patches or another one appearing later is the one becoming unstable first does only dependent indirectly on the minimum value of τ p (x). More important is whether τ p (x) remains low in the near region of the local minimum. The nucleation patch at x/L ≈ 0.75 is in an area of relatively low τ p (x), compared to the other early nucleation patches, which is why it develops faster to the critical size and causes unstable propagation.
This non-local character of the nucleation patches becomes obvious when considering the integral form of Eq. 9 that corresponds to a weighted moving average of τ p (x). In Fig. 4a, we can see that τ n at x/L ≈ 0.75 is considerably lower than at the location of the other early nucleation patches x/L ≈ 0.4 and 0.6. This is why x/L ≈ 0.75 gets critical first and causes unstable slip area growth. Interestingly, x/L ≈ 0.3 is the second most critical point even though the local minimum in τ p is higher than many others in this system. However, τ p (x) remains rather low over an area that approaches h n , and therefore τ n is also low.
In Fig. 4, we compared the prediction of τ cr with measurements from simulations and showed that the prediction works generally well. However, we noticed that for decreasing ξ 0 /h n the discrepancies increase. The theory generally predicts higher τ cr than observed in simulations. We believe that this is caused by nucleation patch interaction, which is neglected in the current theory. The interaction may occur if two nucleation patches are near each other. Nucleation patches cause stress redistribution because the stress inside the patch decreases but, due to equilibrium, it increases in the area near the patch. Therefore, a nucleation patch may cause an increase of the effective stress (compared to the applied load) to another patch, which leads to increased patch size and thus h n is reached already at τ sim cr < τ pred cr . Interestingly, a similar phenomenon has recently been observed in simulations of compressive failure governed by a mesoscopic Mohr-Coulomb criterion, where local damage clusters interact and eventually coalesce to macroscopic failure (Dansereau et al., 2019).
The nucleation patch interaction is likely also the cause for discrepancies observed in the prediction of the nucleation location x cr , as shown in Fig. 4c. While most cases are very well predicted, some simulations present unstable growth that starts from a different location. In these cases, two nucleation patches have very similar critical stress level. However, the (slightly) less critical patch interacts with a neighboring smaller patch and thus becomes unstable at a lower stress level than theoretically expected. This is more likely to occur for systems with low ξ 0 /h n since this increases the likelihood of another local minimum being located close to active nucleation patches. Nevertheless, the critical stress level τ cr remains quantitatively well predicted, as shown Fig. 4b and discussed above, because these secondary effects are minor.
The representative simulation illustrated in Fig. 3 shows that the frictional rupture front, after becoming unstable, does not arrest until it propagated across the entire interface leading to global sliding. This is a general feature of our problem and all our simulations present the same behavior. What is the reason for this run-away propagation? Right after nucleation, the slipping area continues to weaken along its entire length, i.e., δ < d c everywhere. However, after some more growth, it transforms slowly into a frictional rupture front, which is essentially a Griffith's shear crack with a cohesive zone and constant residual strength (Svetlizky and Fineberg, 2014;Svetlizky et al., 2020;Garagash and Germanovich, 2012). The arrest of frictional rupture fronts are governed by an energy-rate balance (Kammer et al., 2015), which states that a rupture continues to propagate as long as the (mode II) static energy release rate G II is larger than the fracture energy, i.e., G II > Γ. In our system, the stress drop ∆τ = τ 0 − τ r is uniform since τ 0 and τ r are uniform. Thus, the static energy release rate grows linearly with rupture length G II ∝ h, and it becomes increasingly difficult to arrest a rupture as it continues to grow. Specifically for our case, we find that G II > Γ max for h/h n 3, where Γ max is Γ from Eq. 14 for τ max p . Hence, once the slipping area reached a size of h/h n 3, nothing can stop it anymore -not even τ max p . For h/h n 3, it is theoretically possible for the slipping area to arrest after some unstable propagation. However, a large increase in τ p (x) would need to occur simultaneously on both side, which is very unlikely, in particular for ξ 0 /h n > 1. Therefore, our assumptions of constant stress drop ∆τ and limited variation of local frictional strength causes arrested rupture fronts to be extremely rare. Nevertheless, since Γ max depends on the probability distribution function f (τ p ), a larger variance, and thus a larger τ max p , would make crack arrest (slightly) more likely.
It is interesting to note, however, that arrest of dynamically propagating slipping areas may occur in other systems. Amon et al. (2017), for instance, showed in their simulations that multiple smaller events nucleate and arrest in order to prepare the interface for a global event. In their system, the initial position along the interface is random as well as the friction properties. Therefore, the available elastic energy, which is the driving force, is random and may fluctuate enough to cause arrest. For the same reason, Geus et al. (2019) observed arrested events of various sizes in simulations with random potential energy along the interface. On the contrary, our system, as outlined above, is characterized by steadily increasing available energy and, thus, behaves differently.
Experimental evidence for arrest of frictional rupture fronts is rather limited. The arrest of confined events observed by Rubinstein et al. (2007) on glassy polymers and by Ke et al. (2018) on granite, is caused by non-uniform loading due to the experimental configuration as demonstrated by Kammer et al. (2015) and Ke et al. (2018Ke et al. ( , 2019. While small scale randomness in the applied shear stress may occur, it does not cause arrest -at best, it may slightly delay or expedite it. Therefore, these experimental observations do not support the presence of any important randomness in the applied shear load; at least at these scales. In much larger systems, such as tectonic plates, randomness in the background stress is likely very important, as discussed in Sec. 5.4.

Interpretation of Monte Carlo Study
In an engineering context, it is usually not enough to know the mean value of a macroscopic property, e.g., the static friction strength, since design criteria are determined based on probability of failure; and risk assessments require failure probability analysis. If the stochastic properties of local interfacial strength τ p (x) are known, the developed theoretical framework in Sec. 4 provides a tool to evaluate the global strength distribution and, hence, the failure probability. However, τ p (x) is not directly observable in experiments (at least so far). In the absence of experimental evidence, the stochastic properties of τ p (x) have been chosen based on physical considerations but our specific parameter choice is arbitrary. This approach is a first approximation that enables us to develop an efficient Monte Carlo method to study the effects of such variations on macroscopic properties. This method may be adapted to more realistic random friction profiles. In the following section we will discuss the choice of each stochastic property of τ p (x) and its effect on the variability of global strength, τ cr .
We assumed that τ p (x) is a random variable following a Beta distribution because it provides simultaneously a non-Gaussian property and well-defined boundaries for minimum and maximum strength, which is physically consistent since mechanical properties are bounded. The parameters of the Beta distribution are chosen such that it is skewed towards the lower bound of τ p . Physically this means that the local interface strength is mostly weak with few strong regions. Under this assumption we observed that the global interface strength τ cr is close the lower bound of τ p and that decreasing correlation length ξ 0 leads to higher τ cr with smaller variation (see Fig. 6). Nucleation is governed by the strength of the weakest region which size equals or exceeds the nucleation length. Hence, if the skewness would be toward high values of τ p -this would corresponds to a mostly strong interface with few weak regions -the variation of τ cr would be larger because of the longer tail at the lower bound. However, the effect of correlation length would remain unchanged: lower ξ 0 would cause higher τ cr with smaller variation.
Further, we assumed that τ p has a power spectral density specified in Eq. 4 with a specific exponent, which affects the memory of the random field. A smaller exponent would result in a flatter decay above the cutoff frequency, and thus generate a field with more high-frequency content. However, when considering equivalent correlation lengths 3 , effects of different assumptions regarding the functional form of the power spectral density are expected to be minor in the probability density of τ cr .
It is noteworthy that there are three relevant length scales: L, h n and ξ 0 . Here, we have not considered the effects of changes in L so far. Based on our theoretical model, we expect that a larger L would result in smaller variance of τ cr . One approach to explore this effect, while avoiding to change the size of the experimental system, could be to modify the normal load, which would affect τ p and τ r , and hence the critical nucleation size h n ∝ (τ p − τ r ) −1 , as shown experimentally by Latour et al. (2013).

Implications for Earthquake Nucleation
In the current study, we are interested in estimating the probability distribution of the macroscopic strength of a frictional interface of given size L. Similar systems but with a focus on other aspects have been studied in order to gain a better understanding of earthquake nucleation. The challenges in studying earthquake nucleation are associated with the size of the system (hundreds of kilometers) and the limited physical access to measure important properties, such as stress state and frictional properties. However, when and how an earthquake nucleates affects directly the average stress drop level ∆τ ≡ τ cr − τ kin , and the earthquake magnitude, which is more easily determined. Hence, there is a need to infer from earthquake magnitude observation back on the fault properties and their variability, to learn about the risk of potential future earthquakes. This is a similar inverse problem as described above.
Previous studies have shown that randomness in simplified models present earthquake magnitudes that follow a power-law distribution (Carlson and Langer, 1989;Ampuero et al., 2006). The simulations presented a large range of magnitudes because the slipping areas arrested, which is the result of randomness in the local stress drop, as discussed in Sec. 5.2. Interestingly, small events were shown to smoothen the stress profile, which reduces the randomness, thus prepared the interface for larger events, as also observed experimentally (Ke et al., 2018). This would suggest that nucleation of larger events tend to be caused by randomness in fault properties rather than (background) stress level, since the stress is getting smoothed. Therefore, our model provides a simple but reasonable tool to study nucleation of medium to large earthquakes.
Our results show that smaller correlation length lead to higher overall strength and more variation. In the context of earthquakes, this would suggest that smaller ξ 0 /h n support larger earthquake magnitudes since nucleation at a higher stress level translates into larger average stress drops, which provides more available energy to release and, thus, makes arrest more difficult. This is complementary to observations by Ampuero et al. (2006) that showed a trend to higher earthquake magnitudes for decreasing standard deviation of the random stress drop field ∆τ(x) while keeping ξ 0 /h n constant.
In addition to the important question on critical stress level for earthquake nucleation, it also remains unclear how the nucleation process takes place. Two possible models (Beroza and Ellsworth, 1996;Noda et al., 2013;McLaskey, 2019) are discussed. The "Cascade" model, where foreshocks trigger each other with increasing size and finally lead to the main earthquake; and the "Pre-Slip" model, which assumes that nucleation is the result of aseismic slow-slip. In our simplistic model, nucleation occurs in a "Pre-Slip" type process, with a long phase of slow slip and an abrupt acceleration after the nucleation patch reached its critical size (see Fig. 3b-right). However, if the amplitude range of our random τ p (x) field was much larger, the likelihood of arrest would increase and, thus, interaction between arrested small events could emerge. This would lead to a nucleation process that resembles more the "Cascade" model. We, therefore, conclude that the type of nucleation process that may occur at a given fault depends on extent of randomness in the local stress and property fields.

Conclusion
We studied the stochastic properties of frictional interfaces considering the nucleation of unstable slip patches. We considered a uniform loading condition and studied the effect of random interface strength, characterized by its probability density and correlation function. Using numerical simulations solving the elastodynamic equations, we demonstrated that macroscopic sliding does not necessarily occur when the weakest point along the interface starts sliding, but when one of possible many slowly slipping nucleation patches reaches a critical length and becomes unstable. We verified that the nucleation criterion originally developed by Uenishi and Rice (2003) predicts well the critical stress leading to global sliding if the criterion is formulated as a minimum of the local strength convolved with the first eigenfunction of the elastic problem. The simulations further showed that increasing correlation lengths of the random interface strength lead to reduced macroscopic static friction. Using the theoretical nucleation criterion, we perform a Monte Carlo study that provided an accurate description of the underlying probability density functions for these observed variations in macroscopic friction. We showed that the probability density function of the global critical strength approaches the probability density of the minimum in the random local strength when the correlation length is much larger than the critical nucleation length. Conversely, a vanishingly small correlation length results in generally higher macroscopic strength with smaller variation. We showed that the presence of precursory dynamic slip events, as in more complex models, is extremely unlikely under the assumption of uniform stress drop. Finally, we discussed discrepancies between the theoretical model and simulations, which suggest that for small correlation lengths the theoretical prediction overestimates the frictional strength, possibly because it neglects interactions between neighboring nucleation patches.