Stochastic Properties of Static Friction

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 $\tau_{0}(t)$ that increases quasi-statically with time (see Fig. 1b). Once $\tau_{0}(t)$ reaches the macroscopic strength of the interface $\tau_{\mathrm{cr}}$, the interface starts to slide and the frictional strength suddenly reduces to its kinetic level $\tau_{\mathrm{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 $\tau_{\mathrm{cr}}$, if divided by the contact pressure, corresponds to the static macroscopic friction coefficient. Similarly, $\tau_{\mathrm{kin}}$ is proportional to the kinetic friction coefficient.

These macroscopic observations depend on the local interface properties, which are the peak strength $\tau_{\mathrm{p}}(x)$, and residual strength $\tau_{\mathrm{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 $\delta(x)$ is local slip, $d_{\mathrm{c}}(x)$ is a characteristic length scale, and $W(x)=(\tau_{\mathrm{p}}(x)-\tau_{\mathrm{r}}(x))/d_{\mathrm{c}}(x)$ is the weakening rate. $H(.)$ is the Heaviside function.

We consider a heterogeneous system with local peak strength $\tau_{\mathrm{p}}(x)$ being a random field, as further described in Sec. 2.2. To reduce complexity of the problem, we assume uniform residual strength¹¹1We use $\partial_{i}$ as short notation for partial derivative with respect to $i$. $\partial_{x}\tau_{\mathrm{r}}=0$ and uniform weakening rate $\partial_{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.

The local peak strength $\tau_{\mathrm{p}}(x)$ is modeled as a stationary non-Gaussian random field with specified cumulative distribution function $F(\tau_{\mathrm{p}})$ and corresponding probability density $f(\tau_{\mathrm{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 $\Phi$ its cumulative distribution, depicted in Fig. 2a-left. $F$ and $\Phi$ are monotonic by definition, so their inverse exist, which can be used to prove that $P\left(\tau_{\mathrm{p}}(x)\leq\tau\right)=F(\tau)$. Further, with $\tau_{\mathrm{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., $\tau_{\mathrm{p}}^{\min}\geq\tau_{\mathrm{r}}$. Second, it maximum value is limited by the material properties. For this reason, we require that $\tau_{\mathrm{p}}\in(\tau_{\mathrm{p}}^{\min},\tau_{\mathrm{p}}^{\max})$, which we achieve by setting $F(\tau_{\mathrm{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}(\xi)$, i.e.,

where $k$ is the angular wave number. We assume that $z(x)$ has a power spectral density

where $\lambda$ is the cutoff frequency, above which the spectral density decays as a power law $\sim k^{-8}$ (see Fig. 2b). The correlation length $\xi_{0}$ is a measure of memory of the random field; the longer $\xi_{0}$ the longer the memory. $\xi_{0}$ is inversely proportional to $\lambda$, and we define²²2The correlation length does not have a precise definition. And alternative definition is $C(\xi_{0})=\exp(-1)$. it as $\xi_{0}=2\pi/\lambda$. Since the correlation function of $z$ is positive, $C_{z}(\xi)>0$, it is not greatly affected by the nonlinear mapping $F^{-1}\circ\Phi$ and $C_{z}(\xi)\approx C_{\tau_{\mathrm{p}}}(\xi)$ (Grigoriu, 1995, p.48) and features such as $\xi_{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.

The samples of random field $\tau_{\mathrm{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\pi j/L$. The fundamental wavelength of the field $2\pi/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 $\tau_{p}(x)$. The modal variance $\sigma_{j}^{2}\propto g(k_{j})$ corresponds to the discrete spectral density, which is normalized to assure that $z$ has unit variance

Due to the discrete representation of $z(x)$, we apply a truncation frequency that is considerably larger than the cutoff frequency $\bar{k}\equiv k_{J}=2.5\lambda$. This ensures that most of the spectral power is preserved:

Further increase in $\bar{k}$ would include additional high frequency modes but with negligible amplitudes. Finally, once $z(x)$ has been generated, we apply the nonlinear mapping $F^{-1}\circ\Phi$ (Eq. 2 and visualized in Fig. 2a) and obtain the random field $\tau_{\mathrm{p}}(x)$. Fig. 1d shows a sample of $\tau_{\mathrm{p}}(x)$ generated using the described procedure with corresponding correlation function and probability density.

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 elasto-dynamic 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~{}\mathrm{GPa}$, Poisson’s ratio of $\nu=0.33$ and density $\rho=1170~{}\mathrm{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 $\tau_{\mathrm{p}}(x)$ as a random field, following the description provided in Sec. 2.2, and constant $\tau_{\mathrm{r}}$. $\tau_{\mathrm{p}}(x)$ follows a Beta distribution with $\alpha=1.5$ and $\beta=3$. We impose a maximum value for relative peak strength of $\max(\tau_{\mathrm{p}}(x)-\tau_{\mathrm{r}})=1.66~{}\mathrm{MPa}$ and minimum value of $\min(\tau_{\mathrm{p}}(x)-\tau_{\mathrm{r}})=0.66~{}\mathrm{MPa}$. Therefore, the random field has a mean value of $\left<\tau_{\mathrm{p}}-\tau_{\mathrm{r}}\right>=1~{}\mathrm{MPa}$ and standard deviation of $0.2~{}\mathrm{MPa}$. We further apply a constant slip-weakening rate of $W=0.5~{}\mathrm{TPa/m}$, which is representative for glassy polymers (Svetlizky et al., 2020). Finally, a slowly increasing uniform stress $\tau_{0}(t)$ is applied along the entire interface.

We use a repetition length of $L=0.1~{}\mathrm{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 $\tau_{\mathrm{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\approx 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 $\tau_{0}(t)=\int_{L}\tau(x,t)\mathrm{d}x$. At the precise moment when the slip patch becomes unstable, marked by a black dot, $\tau_{0}(t)$ starts decreasing rapidly. The maximum value, denoted $\tau_{\mathrm{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.

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 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. 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-\nu)$ for mode II plane-stress ruptures, assuming the stress within $h_{\mathrm{n}}$ has not attained the residual value $\tau_{\mathrm{r}}$ anywhere. Eq. 8 shows that $h_{\mathrm{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 $\tau_{\mathrm{p}}(x)$. Since we have homogeneous elastic solids and a uniform slip-weakening rate $W$, the critical size $h_{\mathrm{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_{\mathrm{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. 21, and can be rewritten in terms of $\tau_{\mathrm{p}}(x)$ and $h_{\mathrm{n}}$ as

where $x$ is the center location of the nucleation patch and $v_{0}(s)\approx(0.925-0.308s^{2})\sqrt{1-s^{2}}$ is the first eigenfunction of the elastic problem. Note that the transformation applied to the argument of $\tau_{\mathrm{p}}(x)$ results in the integral being computed over the critical nucleation patch size $h_{\mathrm{n}}$. Eq. 9 shows that the nucleation stress, which leads to a nucleation patch of size $h_{\mathrm{n}}$, does clearly depend on the shape of $\tau_{\mathrm{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 $\tau_{\mathrm{n}}(x)$ was derived for a contained non-uniformity, for which we know the location. Therefore, $\tau_{\mathrm{n}}(x)$ corresponds to the critical stress of the system. In our system, however, $\tau_{\mathrm{p}}(x)$ is random and multiple nucleation patches might slowly grow. Determining the critical stress $\tau_{\mathrm{cr}}$ of the system requires computing the nucleation stress $\tau_{\mathrm{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 $\tau_{\mathrm{cr}}$ as

For simplicity, we refer to this definition also as $\tau_{\mathrm{cr}}=\min(\tau_{\mathrm{n}}(x))$ and $x_{\mathrm{cr}}=\mathrm{arg}\min(\tau_{\mathrm{n}}(x))$. While it is possible, but not very likely, to have multiple minima of $\tau_{\mathrm{n}}$ with the same amplitude, this does not affect the resulting $\tau_{\mathrm{cr}}$. However, multiple $x_{\mathrm{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.

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 $\tau_{\mathrm{p}}(x)$ generated using the method described in Sec. 2.3. For each of the three different correlation lengths $\xi_{0}/h_{\mathrm{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_{\mathrm{n}}/L=0.034$.

A representative example is shown in Fig. 3. The size of $h_{\mathrm{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 $\tau_{\mathrm{n}}(x)$ (solid blue line) is computed from $\tau_{\mathrm{p}}(x)$ (gray dashed line) using Eq. 9, and $\tau_{\mathrm{cr}}$ is, according to Eq. 10, the minimum of $\tau_{\mathrm{n}}(x)$ (marked by black dot). We find the location of nucleation to be $x_{\mathrm{cr}}/L\approx 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 $\tau_{\mathrm{cr}}^{\mathrm{pred}}$ with the measured value from dynamic simulations $\tau_{\mathrm{cr}}^{\mathrm{sim}}$. We compute $\tau_{\mathrm{cr}}^{\mathrm{pred}}$ as described above with Eq. 10, and as illustrated in Fig. 4a. We further find $\tau_{\mathrm{cr}}^{\mathrm{sim}}=T\partial_{t}\tau_{0}$, where $\partial_{t}\tau_{0}$ is the applied loading rate and $T$ is the time at which $\tau_{0}(t)$ is maximal (see Fig. 3b-left). Comparison of $\tau_{\mathrm{cr}}^{\mathrm{pred}}$ with $\tau_{\mathrm{cr}}^{\mathrm{sim}}$ is shown in Fig. 4b for all $60$ simulations. The results show that the prediction works generally well. For decreasing $\xi_{0}/h_{\mathrm{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 $\tau_{\mathrm{cr}}$ increases with decreasing $\xi_{0}/h_{\mathrm{n}}$.

While the location of nucleation is not relevant for the apparent global strength of our system, we compare the predicted and simulated $x_{\mathrm{cr}}$ for further evaluation of the developed theory. The comparison shown in Fig. 4c uses $x_{\mathrm{cr}}^{\mathrm{pred}}$, as given by Eq. 10 and shown for an example in Fig. 4a, and $x_{\mathrm{cr}}^{\mathrm{sim}}$ 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 $\xi_{0}/h_{\mathrm{n}}=0.25$, the prediction does not work. However, as shown in Fig. 4b, $\tau_{\mathrm{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 $\tau_{\mathrm{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.

The effective frictional strength $\tau_{\mathrm{cr}}=\min(\tau_{\mathrm{n}}(x))$ requires the computation of the nucleation strength $\tau_{\mathrm{n}}(x)$, which involves a convolution of the local peak strength $\tau_{\mathrm{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 $\tau_{\mathrm{p}}(x)$:

where $\tilde{}$ signifies that $\tilde{\tau}_{\mathrm{p}}(x)$ is an approximation of $\tau_{\mathrm{p}}(x)$ and the number of frequencies $J$ is chosen such that the approximation error $|\tilde{\tau}_{\mathrm{p}}-\tau_{\mathrm{p}}|$ is negligible. $\hat{\tau}_{\mathrm{p}}(k_{j})$ is the discrete Fourier transform of $\tau_{\mathrm{p}}(x)$

where $\tau_{\mathrm{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)=\int_{-1}^{+1}e^{-ik_{j}\left(\frac{h_{\mathrm{n}}}{2}s+x\right)}v_{0}(s)\mathrm{d}s$ is the modal convolution term, which, being independent of the sample specific functional form of $\tau_{\mathrm{p}}(x)$, can be pre-computed. This formulation allows for efficient and precise evaluation of the effective frictional strength $\tau_{\mathrm{cr}}=\min\tau_{\mathrm{n}}(x)$ for a large number of samples $N=10,000$, such that the probability distribution $f(\tau_{\mathrm{cr}})$ and its evolution as function of the correlation length $\xi_{0}$ can be accurately studied.

Prior to presenting the numerical results we provide some intuition of the effect of correlation length on the nucleation strength $\tau_{\mathrm{n}}$ based on probabilistic arguments. By exploiting the stationarity of $\tau_{\mathrm{p}}$ and $\tau_{\mathrm{n}}$ it is possible to derive an analytical expression of the expectation of the nucleation strength $\mathrm{E}[\tau_{\mathrm{n}}]$ and its variance $\mathrm{Var}[\tau_{\mathrm{n}}]$ as function of the corresponding statistical properties of local strength, $\mathrm{E}[\tau_{\mathrm{p}}]$, $\mathrm{Var}[\tau_{\mathrm{p}}]$ and $\xi_{0}/h_{\mathrm{n}}$ (see B).

One interesting finding is that the expectation is not affected by $\xi_{0}/h_{\mathrm{n}}$: $\mathrm{E}[\tau_{\mathrm{n}}]=\mathrm{E}[\tau_{\mathrm{p}}]$ (see derivation in Eq. 23). The expression for $\mathrm{Var}[\tau_{\mathrm{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. 24). For perfect correlation, i.e., $\xi_{0}/h_{\mathrm{n}}=\infty$, $C(.)$ becomes a constant, thus $\mathrm{Var}[\tau_{\mathrm{n}}]=\mathrm{Var}[\tau_{\mathrm{p}}]$. Additionally, in the limit of $\xi_{0}\ll h_{\mathrm{n}}$, the double integral in Eq. 24 scales with $\xi_{0}/h_{\mathrm{n}}$, thus $\mathrm{Var}[\tau_{\mathrm{n}}]\propto\mathrm{Var}[\tau_{\mathrm{p}}]{\xi_{0}}/{h_{\mathrm{n}}}$ (see derivation in Eq. 26).

We consider a range of correlation lengths $\xi_{0}/h_{\mathrm{n}}=\{0.25,0.5,1.0,2.0\}$, while all other properties remain constant. Fig. 5a-left shows one sample of $\tau_{\mathrm{p}}(x)$ for each considered $\xi_{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 $\sigma_{j}^{2}$, corresponding to the different $\xi_{0}$. For this reason, all shown samples have a similar spatial evolution and the effect of varying $\xi_{0}$ can be clearly observed. By definition, all $\tau_{\mathrm{p}}(x)$ samples are drawn from the same probability distribution $f(\tau_{\mathrm{p}})$ (see Fig. 5a-center). Decreasing $\xi_{0}$, moves the probability density of its minimum $f(\min(\tau_{\mathrm{p}}))$ towards the lower bound $\tau_{\mathrm{p}}^{\mathrm{min}}=0.66\langle\tau_{\mathrm{p}}\rangle$ (see Fig. 5a-right), because with lower correlation lengths it is more likely to visit a broad range of $\tau_{\mathrm{p}}$ values.

Fig. 5b-left shows the corresponding nucleation strength $\tau_{\mathrm{n}}(x)$ for each of the local frictional strength fields $\tau_{\mathrm{p}}(x)$ presented in Fig. 5a, computed using Eq. 13. As mentioned before, $\tau_{\mathrm{n}}$ is essentially a weighted moving average of $\tau_{\mathrm{p}}$ with window size $h_{\mathrm{n}}$ (see Eq. 9). Thus, most of the high frequency content of $\tau_{\mathrm{p}}$ disappears and the effect of $\xi_{0}$ on $\tau_{\mathrm{n}}$ is more subtle. One interesting feature is in the minima and maxima of $\tau_{\mathrm{n}}$: increasing $\xi_{0}$ causes lower minima and higher maxima, because the moving average is effectively computed over an approximately constant field $\tau_{\mathrm{n}}\approx\tau_{\mathrm{p}}$. Inversely, decreasing $\xi_{0}$ causes the opposite effect and $\tau_{\mathrm{n}}\approx\langle\tau_{\mathrm{p}}\rangle$.

This effect is more clearly visible by considering the distribution $f(\tau_{\mathrm{n}})$ shown in Fig. 5b-right. Increasing $\xi_{0}$ effectively puts more weight onto the tails of $f(\tau_{\mathrm{n}})$ (see $\xi_{0}/h_{\mathrm{n}}=2.0$ in Fig. 5b-right), and in the limiting case of $\xi_{0}/h_{\mathrm{n}}\rightarrow\infty$ the distribution of $\tau_{\mathrm{n}}$ will be the same as the one of $\tau_{\mathrm{p}}$ (analogous to Eq. 25). On the other hand, decreasing $\xi_{0}$ puts weight on its mean $\langle\tau_{\mathrm{p}}\rangle$, making $f(\tau_{\mathrm{n}})$ similar to a Gaussian (see $\xi_{0}/h_{\mathrm{n}}=0.25$ in Fig. 5b-right) with variance proportional to $\xi_{0}$ (see Eq. 26). In the limit $\xi_{0}/h_{\mathrm{n}}\rightarrow 0$ the distribution of $\tau_{\mathrm{n}}$ becomes a Dirac-$\delta$ centered at $\langle\tau_{\mathrm{p}}\rangle$. The described dependence of $f(\tau_{\mathrm{n}})$ on $\xi_{0}$ confirms the previously stated statistical arguments (see B for derivation).

Because $f(\tau_{\mathrm{p}})$ is skewed towards the lower bound of $\tau_{\mathrm{p}}$ so is $f(\tau_{\mathrm{n}})$; the larger $\xi_{0}$ the larger the skewness. For $\tau_{\mathrm{cr}}$ this effect is amplified by the fact that $\tau_{\mathrm{cr}}=\min(\tau_{\mathrm{n}}(x))$ as depicted in Fig. 5c, causing $\langle\tau_{\mathrm{cr}}\rangle$ to decrease with increasing $\xi_{0}$. As noted in Sec. 3.3, the location where the critical instability occurs $x_{\mathrm{cr}}$ is uniformly distributed over the entire domain as shown in Fig. 5d and is independent on $\xi_{0}$.

We further analyze the effect of $\xi_{0}$ on the probability distribution of $\tau_{\mathrm{n}}$ and $\tau_{\mathrm{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_{\xi_{0}/h_{\mathrm{n}}\rightarrow 0}\tau_{\mathrm{n}}(x)=\langle\tau_{\mathrm{p}}\rangle$, because the moving average in computing $\tau_{\mathrm{n}}(x)$ is evaluated over a window $h_{\mathrm{n}}$ that appears infinite compared to $\xi_{0}$ (see Fig. 6a). Consequently, the effective strength also tends towards the mean of the local strength: $\lim_{\xi_{0}/h_{\mathrm{n}}\rightarrow 0}\tau_{\mathrm{cr}}=\langle\tau_{\mathrm{p}}\rangle$ (see Fig. 6b). Conversely, if ${\xi_{0}\gg h_{\mathrm{n}}}$, the moving average is computed over a window $h_{\mathrm{n}}$ which vanishes, and thus (9) becomes the identity: $\lim_{\xi_{0}/h_{\mathrm{n}}\rightarrow\infty}\tau_{\mathrm{n}}(x)=\tau_{\mathrm{p}}(x)$. In this case, the effective strength will be more likely to be close to the actual lower bound of the distribution $\lim_{\xi_{0}/h_{\mathrm{n}}\rightarrow\infty}\tau_{\mathrm{cr}}=\min(\tau_{\mathrm{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 $\xi_{0}/h_{\mathrm{n}}$ by reporting the mean effective friction strength (see inset in Fig. 6b). For $\xi_{0}/h_{\mathrm{n}}\geq 0.5$ simulations and theory coincide. However, for $\xi_{0}/h_{\mathrm{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 $\xi_{0}/h_{\mathrm{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.

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., $\tau_{\mathrm{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 $\sim 10-100~{}\mathrm{mm}$ (Ben-David and Fineberg, 2011; Latour et al., 2013) and surface roughness lengths of $\sim 1~{}\mu\mathrm{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 $\delta$, slip rate $\partial_{t}\delta$, 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_{\mathrm{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 rate-weakening 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 $\sim 1~{}\mathrm{MPa}$ for glassy polymers at normal pressure of $\sim 5~{}\mathrm{MPa}$ (Svetlizky et al., 2020) and, similarly, $\sim 1~{}\mathrm{MPa}$ weakening for granite at normal pressure of $\sim 6~{}\mathrm{MPa}$ (Kammer and McLaskey, 2019), we considered a model system with strong frictional weakening. However, we neglect the complexity of rate-and-state friction, as extensively demonstrated by Ray and 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 $\Gamma$, which is the energy dissipated by the weakening process, i.e., the triangular area $(\tau_{\mathrm{p}}-\tau_{\mathrm{r}})d_{\mathrm{c}}/2$ in Fig. 1c:

Since $W$ is constant in our system $\Gamma$ varies with $(\tau_{\mathrm{p}}(x)-\tau_{\mathrm{r}})^{2}$. A correlation between $\Gamma(x)$ and $\tau_{\mathrm{p}}(x)$ can be expected given that any point that is stronger, i.e., increased $\tau_{\mathrm{p}}$, is likely to dissipate more energy as well, i.e., increased $\Gamma$.

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 $\tau_{\mathrm{p}}$, $\tau_{\mathrm{r}}$ and $d_{\mathrm{c}}$.

Finally, we note that by assuming a periodic system, we neglect possible boundary effects. We expect that the boundary would locally reduce $\tau_{\mathrm{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 $\tau_{\mathrm{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_{\mathrm{n}}$ and, therefore, $f(\tau_{\mathrm{cr}})$ will tend towards the periodic solution for $h_{\mathrm{n}}/L\rightarrow 0$. Verification would require a large number of numerical simulations, which is beyond the scope of this work.

The simulations have shown that uniform $\tau_{0}$ and random $\tau_{\mathrm{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 $\tau_{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 $\tau_{\mathrm{p}}(x)$ profile shown in Fig. 3c presents three local minima with approximately the same value, i.e., at $x/L\approx 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 $\tau_{\mathrm{p}}(x)$. More important is whether $\tau_{\mathrm{p}}(x)$ remains low in the near region of the local minimum. The nucleation patch at $x/L\approx 0.75$ is in an area of relatively low $\tau_{\mathrm{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 $\tau_{\mathrm{p}}(x)$. In Fig. 4a, we can see that $\tau_{\mathrm{n}}$ at $x/L\approx 0.75$ is considerably lower than at the location of the other early nucleation patches $x/L\approx 0.4$ and 0.6. This is why $x/L\approx 0.75$ gets critical first and causes unstable slip area growth. Interestingly, $x/L\approx 0.3$ is the second most critical point even though the local minimum in $\tau_{\mathrm{p}}$ is higher than many others in this system. However, $\tau_{\mathrm{p}}(x)$ remains rather low over an area that approaches $h_{\mathrm{n}}$, and therefore $\tau_{\mathrm{n}}$ is also low.

In Fig. 4, we compared the prediction of $\tau_{\mathrm{cr}}$ with measurements from simulations and showed that the prediction works generally well. However, we noticed that for decreasing $\xi_{0}/h_{\mathrm{n}}$ the discrepancies increase. The theory generally predicts higher $\tau_{\mathrm{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_{\mathrm{n}}$ is reached already at $\tau_{\mathrm{cr}}^{\mathrm{sim}}<\tau_{\mathrm{cr}}^{\mathrm{pred}}$. 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_{\mathrm{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 $\xi_{0}/h_{\mathrm{n}}$ since this increases the likelihood of another local minimum being located close to active nucleation patches. Nevertheless, the critical stress level $\tau_{\mathrm{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., $\delta<d_{\mathrm{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_{\mathrm{II}}$ is larger than the fracture energy, i.e., $G_{\mathrm{II}}>\Gamma$. In our system, the stress drop $\Delta\tau=\tau_{0}-\tau_{\mathrm{r}}$ is uniform since $\tau_{0}$ and $\tau_{\mathrm{r}}$ are uniform. Thus, the static energy release rate grows linearly with rupture length $G_{\mathrm{II}}\propto h$, and it becomes increasingly difficult to arrest a rupture as it continues to grow. Specifically for our case, we find that $G_{\mathrm{II}}>\Gamma_{\mathrm{max}}$ for $h/h_{\mathrm{n}}\gtrapprox 3$, where $\Gamma_{\mathrm{max}}$ is $\Gamma$ from Eq. 14 for $\tau_{\mathrm{p}}^{\mathrm{max}}$. Hence, once the slipping area reached a size of $h/h_{\mathrm{n}}\gtrapprox 3$, nothing can stop it anymore – not even $\tau_{\mathrm{p}}^{\mathrm{max}}$. For $h/h_{\mathrm{n}}\lessapprox 3$, it is theoretically possible for the slipping area to arrest after some unstable propagation. However, a large increase in $\tau_{\mathrm{p}}(x)$ would need to occur simultaneously on both side, which is very unlikely, in particular for $\xi_{0}/h_{\mathrm{n}}>1$. Therefore, our assumptions of constant stress drop $\Delta\tau$ and limited variation of local frictional strength causes arrested rupture fronts to be extremely rare. Nevertheless, since $\Gamma_{\mathrm{max}}$ depends on the probability distribution function $f(\tau_{\mathrm{p}})$, a larger variance, and thus a larger $\tau_{\mathrm{p}}^{\mathrm{max}}$, 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. (2018, 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.

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 $\tau_{\mathrm{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, $\tau_{\mathrm{p}}(x)$ is not directly observable in experiments (at least so far). In the absence of experimental evidence, the stochastic properties of $\tau_{\mathrm{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 $\tau_{\mathrm{p}}(x)$ and its effect on the variability of global strength, $\tau_{\mathrm{cr}}$.

We assumed that $\tau_{\mathrm{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 $\tau_{\mathrm{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 $\tau_{\mathrm{cr}}$ is close the lower bound of $\tau_{\mathrm{p}}$ and that decreasing correlation length $\xi_{0}$ leads to higher $\tau_{\mathrm{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 $\tau_{\mathrm{p}}$ – this would corresponds to a mostly strong interface with few weak regions – the variation of $\tau_{\mathrm{cr}}$ would be larger because of the longer tail at the lower bound. However, the effect of correlation length would remain unchanged: lower $\xi_{0}$ would cause higher $\tau_{\mathrm{cr}}$ with smaller variation.

Further, we assumed that $\tau_{\mathrm{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$C(\xi_{0})=e^{-1}$, effects of different assumptions regarding the functional form of the power spectral density are expected to be minor in the probability density of $\tau_{\mathrm{cr}}$.

It is noteworthy that there are three relevant length scales: $L$, $h_{\mathrm{n}}$ and $\xi_{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 $\tau_{\mathrm{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 $\tau_{\mathrm{p}}$ and $\tau_{\mathrm{r}}$, and hence the critical nucleation size $h_{\mathrm{n}}\propto(\tau_{\mathrm{p}}-\tau_{\mathrm{r}})^{-1}$, as shown experimentally by Latour et al. (2013).

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 $\langle\Delta\tau\rangle\equiv\tau_{\mathrm{cr}}-\tau_{\mathrm{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 $\xi_{0}/h_{\mathrm{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 $\Delta\tau(x)$ while keeping $\xi_{0}/h_{\mathrm{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 $\tau_{\mathrm{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.