On the speed of fast and slow rupture fronts along frictional interfaces

The transition from stick to slip at a dry frictional interface occurs through the breaking of the junctions between the two contacting surfaces. Typically, interactions between the junctions through the bulk lead to rupture fronts propagating from weak and/or highly stressed regions, whose junctions break first. Experiments find rupture fronts ranging from quasi-static fronts with speeds proportional to external loading rates, via fronts much slower than the Rayleigh wave speed, and fronts that propagate near the Rayleigh wave speed, to fronts that travel faster than the shear wave speed. The mechanisms behind and selection between these fronts are still imperfectly understood. Here we perform simulations in an elastic 2D spring--block model where the frictional interaction between each interfacial block and the substrate arises from a set of junctions modeled explicitly. We find that a proportionality between material slip speed and rupture front speed, previously reported for slow fronts, actually holds across the full range of front speeds we observe. We revisit a mechanism for slow slip in the model and demonstrate that fast slip and fast fronts have a different, inertial origin. We highlight the long transients in front speed even in homogeneous interfaces, and we study how both the local shear to normal stress ratio and the local strength are involved in the selection of front type and front speed. Lastly, we introduce an experimentally accessible integrated measure of block slip history, the Gini coefficient, and demonstrate that in the model it is a good predictor of the history-dependent local static friction coefficient of the interface. These results will contribute both to building a physically-based classification of the various types of fronts and to identifying the important mechanisms involved in the selection of their propagation speed.


I. INTRODUCTION
The onset of sliding at a frictional interface occurs through the breaking of the many contacts that were preventing the relative motion of the surfaces. When a single contact breaks, the stress it bore is redistributed to its neighbors, bringing them closer to or past their load-bearing capacity. In extended frictional interfaces (i.e. larger than a characteristic correlation length scale), this process can lead to propagating ruptures -rupture fronts. The recent direct observation of rupture fronts in laboratory friction experiments (see e.g. [1][2][3][4][5]) have provided new insights and opened new questions. It was found, for instance, that not all fronts span the entire interface [6,7]. The selection of the propagation length of the fronts has been intensely investigated [8][9][10][11][12][13][14][15]. It was also found that the fronts can propagate at a variety of speeds, either quasi-statically [4,5,16], at speeds close to that of surface waves (sub-Rayleigh) [3,17,18], at speeds faster than the shear wave speed c s (supershear) [1,3,[17][18][19], or at speeds one or two orders of magnitude slower than the Rayleigh speed (slow) [3,17,20,21]. In this context, the present paper is mainly devoted to the study of the mechanisms responsible for front speed selection.
Experiments to shed light on the nature of the rupture fronts have been performed. The authors of [17,22] * j.k.tromborg@fys.uio.no † julien.scheibert@ec-lyon.fr placed arrays of sensors close to the interface and used continuum theory to infer the properties of the elastic fields at the interface. Svetlizky and Fineberg [23] recently showed that the dynamic fields associated with sub-Rayleigh fronts are consistent with the ones predicted by linear elastic fracture mechanics [24]. The authors of [5,21,25] used microtextured surfaces where each contact can be tracked individually to follow the rupture fronts directly at the individual micro-contact level. Common to all these experiments is that they seek to improve measurements of the behavior at the very interface by increasing the spatial and temporal resolution. These experiments can be usefully complemented by computer simulations, which can provide, for example, simultaneous access to shear stress, normal stress, local contact strength and front propagation, a combination which remains hard to access experimentally.
To represent the propagation of the front separating a stuck part of the interface from a slipping one, models of the transition to sliding need to include at least one level of discretization of the macroscopic interface. Depending on the model, the stress transfer at this so-called mesoscopic scale may be treated with a 1D [7,8,10,12,[26][27][28][29][30][31][32], 2D [9,11,[33][34][35][36] or fully 3D discretization. The simplest bulk model that includes spatial and temporal structure in the transition to sliding is probably the 1D springblock model, which has been popular in the earthquake literature (see e.g. [37,38]) since the spring-block experiment and simulation by Burridge and Knopoff [39]. As spatiotemporal data for the onset of sliding became available in laboratory friction experiments, the 1D spring-block model was also applied in the friction literature [7,10,12,26,29,30,40]. The main limitations of 1D models are their inability to accurately reproduce the stress fields that arise in the experiments and the lack of a physical length scale [41,42] unless such a length scale is introduced in the friction law (e.g. [10]).
To better reproduce the experimental loading conditions and how they translate to heterogeneous shear and normal stress fields at the frictional interface, we used, in [9], a 2D spring-block model, which can be shown to reproduce 2D linear elasticity [43]. With Amontons-Coulomb friction this model agreed well with experiments for static measures related to the onset of sliding, such as for the length of precursors and the evolution of the normal stress along the interface, but the model did not capture the full dynamics of the rupture fronts. Radiguet et al. [11] studied the memory of the stress state through the passage of multiple ruptures in a visco-elastic 2D finite element model with a slip-weakening friction law. They too focused on successive precursors rather than the dynamics of each rupture event. Kammer et al. [34] studied the properties of fast rupture front propagation in a 2D finite element model with a static + velocity weakening dynamic friction law. Matsukawa and coworkers [35,44] studied how the normal force and the size of the interface influence the effective macroscopic static friction coefficient, using a 2D finite element model with a local velocity-weakening Amontons-Coulomb friction law. While each of the above-mentioned models was able to reproduce some aspects of the properties of rupture fronts, collectively, they do not provide a single framework able to account for the full richness of front dynamics that was observed experimentally. In particular, the self-selection of front type leading to sub-Rayleigh, slow and supershear fronts, as well as the transitions between them, was missing.
One common feature of these 2D models was that they used continuum laws to model friction. Such laws, from the Amontons-Coulomb description, through slipweakening [45][46][47] and velocity-weakening [38], to rateand-state friction laws [48][49][50][51][52], reproduce the robust average behavior of the myriad micro-junctions that make up each mesoscopic region of the frictional interface. However, by their nature, they do not explain how the individual micro-junctions evolve and interact to produce the overall friction behavior. To approach this question, numerical [26,29,36,[53][54][55] and analytical [56][57][58][59] models have been made that explicitly include a set of junctions, each representing one or a few microscopic contacts.
In [36] we combined an asperity model of the friction at the interface with a 2D elastic solver in order to accurately reproduce the experimental loading conditions used in [3]. In the asperity model we included a time scale inspired by the time scale identified in the same experimental system [60] that in the model controls the healing of the interface back to a fully pinned state after slipping. This combined model produced spatiotemporal features of the rupture dynamics very similar to those observed in the experiment. In particular, we reproduced the abrupt transitions between fast and slow front propagation, which can be understood from the underlying fast and slow slip mechanisms. Here, we use the same model to gain significant additional insights into the relation between the microscale junction dynamics, the mesoscale slow and fast slip dynamics, and the macroscale friction dynamics and front propagation. The paper is structured as follows. In Section II we introduce the model and highlight some similarities and differences with other popular friction models. In Section III we introduce the two principal speeds in the problem -the slip and front speeds.
In Section IV we identify a signature of slow fronts in the macroscopic loading curve, and demonstrate the inertial nature of fast slip and fast fronts. We map out the initial conditions that lead to fast and to slow fronts. In Section V we investigate speed selection in the model. We consider successively the spatial extent of front speed transients, the influence of local stress and strength on front speed within each type of front (fast or slow), and the relationship between slip speed and front speed. In Section VI we describe how a rupture event sets up the state of the interface, which in turn determines the properties of the next rupture event. Section VII is a brief discussion.

A. Bulk modeling
We consider the rough frictional interface between a horizontal track and a thin linearly elastic slider (Fig. 1a). The slider has mass M and sizes L and H in the horizontal (x) and vertical (z) directions, respectively. We present the values of all parameters in Table I in Appendix A. The bulk elastodynamics of the slider are solved using a square lattice of N = N x N z point masses (of mass m = M/N ) connected by internal springs (Fig. 1b) [9,43]. Rotational degrees of freedom within the bulk are not included, but we will still designate the point masses as blocks, as is common in models of this type. Blocks are coupled to their four nearest neigh-bours and their four next-nearest neighbors by springs of equilibrium lengths l = L/(N x − 1) = H/(N z − 1) and √ 2l and stiffnesses k and k/2, respectively, giving an isotropic elastic model with Poisson's ratio 1/3. The spring force exerted on block i by block j is thus k ij (r ij − l ij ) ∆xij rij when blocks are connected, 0 otherwise, where x = (x, z), ∆x ij = x j − x i , r ij = |∆x ij | and k ij and l ij are the stiffness and equilibrium length of the spring connecting blocks i and j. Block oscillations are damped by introducing a viscous force acting against the relative motion of connected blocks. For every block i the viscous force from block j is η(ẋ j −ẋ i ) if the blocks are connected, zero otherwise. We chose the coefficient η = √ 0.1km so that blocks are underdamped and event-triggered oscillations die out well before the next event. The top blocks are submitted to uniformly distributed, time-independent vertical forces F N Nx . The bottom blocks lie on an elastic foundation of modulus k f = k/2, i.e. each block is submitted to a vertical force of amplitude p i = k f |z i | if z i < 0 or 0 otherwise, where z i is the vertical displacement of block i. Both vertical boundaries are free, except for a horizontal driving force F T = K(V t − x h ) applied on the left-side block situated at height h above the interface, where x h is the xdisplacement of this block. This models a pushing device of stiffness K driven at a small constant velocity V . The 2N equations of motion are solved simultaneously using a leapfrog / velocity Verlet integrator [61] on a uniform temporal grid of resolution ∆t.

B. Interface modeling
The net contact between two solids generically consists of a large number of stress bearing micro-junctions. The properties of these junctions depend on the type of interface. For rough solids, each micro-junction corresponds to a micro-contact between asperities on the opposing surfaces, whereas for smoother surfaces the junctions can be solidified patches of an adsorbate layer [55]. We include in our model the following three physical aspects of the junction behavior. 1) A micro-junction in its pinned state behaves elastically and can bear a shear force f T , provided f T remains smaller than a threshold f thres . When f thres is reached, a local fracture-like event occurs, and the junction enters a slipping state. 2) In the slipping state, the micro-junction moves with the slider's surface. The physical picture can be e.g. the microslipping of micro-asperities in contact or the fluidization of an adsorbate layer. During slip, the micro-junction sustains some residual force f T = f slip , with f slip smaller than f thres . 3) Slipping micro-junctions have a certain probability to disappear or relax. For example, a microcontact disappears when an asperity moves away from its antagonist asperity by a typical distance equal to the mean size of micro-contacts, as classically considered for slow frictional sliding, e.g. in rate-and-state friction laws. However, another picture may arise from the sudden re-lease of energy when pinned junctions break. This energy will dissipate partly as radiated elastic waves, but mainly as heat in the region around the micro-junction [60]. The rise in temperature will significantly increase the rate of a thermally activated relaxation of the slipping microjunction Such thermal processes will classically lead to time-controlled shear-relaxations, during the time interval needed for the interface to cool down. Typical values of this cooling time, which correlates to the slip dynamics of the interface just after a rupture event [60], have recently been estimated to be in the sub-millisecond range [62], i.e. of the same order as the duration of rupture events themselves. We thus expect this short time scale to be highly relevant for rupture front dynamics and include it as a crucial ingredient of our model.
The physical aspects described above have been modeled in a simple way using the following assumptions. The multi-contact nature of the interface is modeled through an array of N s tangential springs representing individual micro-junctions, attached in parallel to each interfacial block (Fig. 1c) [26,29]. The individual spring behavior is the same as in [36] (see SI Methods therein for additional details). A spring pinned to the track stretches linearly elastically as the block moves, acting with a tangential force f T on the block. When the force reaches the static friction threshold f thres (we neglect aging, so that f thres is time independent), the micro-junction ruptures and the spring becomes a slipping ideal plastic spring acting with a dynamic friction force f T = f slip . We take f thres and f slip proportional to the normal force p on the corresponding block. This ascribes the pressure dependence of friction to the individual forces bore by a constant number N s of springs rather than to a pressure dependent number of springs per block. After a random time t R drawn from a distribution T (t R ), the slipping spring relaxes. It is replaced immediately by a pinned, unloaded spring (f new = 0) representing a new junction formed elsewhere and a new cycle starts. Here we use T (t R ) as a simplified way of modeling the distribution of times after which micro-junctions relax. Due to the variety and the complexity of the underlying thermal processes, we did not try to derive T (t R ) for a specific situation. Rather, we chose to model T (t R ) in the simplest way, as a Gaussian with average time t R and width δt R , modified so that negative times are excluded [63]. The shape of T (t R ) is not crucial: we obtain qualitatively similar results with an exponential distribution. The width of T (t R ) is the only source of randomness in our model and causes the interface springs of a block to evolve differently from each other.

C. Global behavior
For our chosen set of parameters, when loaded from the side the system enters regular stick-slip, with alternating events showing partial and full breaking of the interface (small and large drops in F T , see Fig. 2 two types of events will later be denoted as partial slip and full/global sliding events, respectively. In the model, as in experiments, the behavior on larger scales emerges from the interactions on the scales below. Because the junctions are not all at their force threshold at the same time, the effective static friction of each mesoblock is less than the sum of the individual junction thresholds. Similarly, the effective static friction of the entire slider is less than the sum of the mesoblock static friction levels. As seen in Fig. 2, the global static friction during stick-slip is max(F T /F N ) ≈ 0.19.

D. Relationship to other models
A spring-block (sometimes called mass-spring) discretization of the bulk elasticity is particularly convenient for models where the friction is described as an ensemble of micro-junctions rather than a continuum law, because the blocks provide natural units on which to couple the frictional and the bulk elastic behavior. Like finite element (FEM) and finite difference (FDM) methods, spring-block discretization satisfies the equations of linear elasticity. In particular, longitudinal (P) and shear (S) waves in the bulk propagate with the correct speeds and the right reflection and refraction properties [43]. We checked that our model reproduces the expected bulk wave speeds.
The choice of the spatial resolution is made according to the following physical arguments. First, as discussed by e.g. Persson [56], Caroli and Nozières [64] and Braun et al. [65], below a characteristic length scale λ, the so-called elastic screening length, the interface behaves rigidly. λ is thus the maximum lattice constant allowing for a correct representation of the elasticity of the interface. Second, the frictional behavior of each block is then a statistical average over the many micro-junctions connected to the interface region it represents. This statistical approach is increasingly relevant for coarser meshes where each block involves more junctions. One thus looks for the largest possible lattice constant. Combining both requirements, λ appears as the natural lattice constant for such spring-block models. For a linearly elastic rough interface, λ ∼ d 2 /a, with a the typical lateral size of microcontacts and d the typical distance between them. For micrometer-ranged roughnesses, we expect a ∼ 1 µm and d ∼ 10 − 100 µm, yielding λ ∼ 0.1 − 10 mm.
The rate-and-state picture of friction, which includes displacement-controlled disappearance of micro-contacts, has proved to be adequate for slow (typically up to 100 µm/s range) sliding in a variety of materials. Such slow velocities imply a negligible temperature rise of the interface and thus a slipping state the duration of which is of purely geometrical origin, i.e. it is controlled by a length scale of the order of the mean micro-contact size.
Here, we focus on a drastically different situation, in which the transition from static to kinetic friction is extremely short (millisecond range) and is accompanied by fast slip (100 mm/s range). As recognized in recent experiments by [60,62,66], the sudden rupture of the interface and its subsequent slip will generate a significant heating of the interface, sufficient to melt the broken micro-asperities. In these severe conditions, precise knowledge about the micro-contact behavior at the millisecond time scale is currently lacking.
We emphasize that experimental data did show that the dynamics at the onset of sliding involves a transition from fast to slow slip which occurs after a constant time, rather than a constant displacement [see 60]. It is therefore natural to propose an alternative picture that incorporates the possibility that transitions between the slipping and the pinned states can be controlled by a time. Note that this timescale controls the fast dynamics with which the interface comes back to a fully pinned state after slipping. It thus drastically differs from the classical, longer time scale for aging, also found in [60], and which controls the slow strength recovery of the interface when it has already come back to rest. A complete model including in the junction law a length scale for the removal of the junctions, a time scale for aging [67,68], and the time scale we study here, is beyond the scope of this work.
In terms of modeling approach, let us also stress the fact that a number of reference models from the literature considered, before us, time-controlled transitions between micro-junction states [see e.g. 26,29,55,59,69].
As shown in Thøgersen et al. [59], our friction law is actually one particular case of a more general family of models. Note that this family includes, as another particular case, the purely slip-dependent friction laws that arise when junctions reform immediately after breaking, as studied in [57,58].

III. SLIP AND FRONT DYNAMICS
In the transition from the stuck to the slipping interface there are two types of speed to consider: the material slip speed and the rupture front speed. We will later show that they are closely related, but we first define both in this section, which serves as necessary background information for the new results that follow in Sections IV and V.

A. Rupture front characteristics
With the driving force applied on the trailing edge of the slider, the blocks near the trailing edge are the first to reach their effective static friction thresholds. A block that slips increases the load on its neighbors, which can start to slip in turn. The rupture front tip, i.e. the boundary between a region of stuck blocks and a region of slipping blocks, then propagates away from the nucleation point.
In the interaction law the distinction between pinned and slipping states is made on the junction rather than the block level. We therefore define slipping on the block level to mean that a certain fraction of the block's junctions are in the slipping state. This criterion is robust to the choice of threshold fraction, because as is seen in Fig. 3 and other figures (8,19), it is typical for a block to go from having nearly all its junctions pinned to having nearly all of them broken in a time short compared to the other time scales in the simulation. For the event in Fig. 3 the time to go from 80% of springs pinned to 20% of springs pinned is approximately 0.03 ms in the fast part of the front and 0.3 ms in the slow part of the front. We have used a threshold value of 30%, that is, the start time of block slipping is taken as the instant when the fraction of pinned junctions dropped below 30%. Figure 3a shows a rupture front whose speed v c changes from fast (v c ∼ c s /3) to slow (v c ∼ c s /100) and back to fast again. The left-travelling front that starts when the primary front is reflected from the leading edge re-breaks the junctions that had healed behind the front tip. By defining the transition to block sliding as above, the location of the front tip in time can be measured. The local front speed is then determined from the time it takes the front to travel a fixed number of lattice constants. Because of disorder in the junction state along the interface remaining from earlier events, the propagation time for the front to move from one block to the next can vary significantly. We have found that using only the end-points in a 5 blocks wide moving stencil (Appendix B) gives the best balance between robustness and spatial resolution in the calculated front speed. The results are shown in Fig. 3b.
B. Block slip dynamics and a slow slip mechanism Figure 4a shows the slip dynamics of a block in a partial slip event and a block in a full sliding event. We find both fast and slow slip regimes of the motion, in excellent agreement with the experiments reported in [60]. The initial fast slip regime begins after the passage of a fast rupture front. In both these cases the fast slip is followed by the block coming nearly to rest (no visible increase in net slip between 450 and 550 µs) and then by a slow slip regime with roughly linear increase of slip vs time, i.e. constant slip speed. The slow slip regime can end in two ways. In full sliding events, slow slip changes back to fast slip when the slider enters the full sliding regime. For arresting events, the slow slip regime ends when the block comes to rest.
The initial fast slip regime corresponds to an inertial motion of the mesoblock when a large number of junctions break in a short time interval as the rupture front passes by. The net friction force is rapidly reduced, bringing the block out of mechanical equilibrium. The result    is a large positive acceleration (in the direction of the net force due to the neighboring blocks). The inertial nature of this motion is demonstrated in Fig. 6.
The subsequent slow slip observed in the model has a different physical origin that was explained in [36,59]. This slow slip mechanism is illustrated in Fig. 4b. When f new < f slip the net friction force on a block is reduced slightly whenever a spring leaves the slipping state, yielding a small positive acceleration as the net friction drops below the net external force from neighboring blocks. The friction reduction is soon balanced by the changes in the forces from the pinned junctions and the external forces on the block as it slowly moves. This slow slip mechanism is present as long as some junctions are going from the slipping to the pinned state and f new < f slip , but it is masked by the fast slip while the fast slip lasts. In [36], we discussed in detail the dependence of slow slip speed on model parameters. We also showed that, in the model, slow slip is a sufficient condition for slow fronts to be possible, irrespective of the precise physical origin of slow slip on the block level (see the two alternative models tested in [36]).

IV. FRONT TYPE RESULTS
In this section we discuss the conditions under which we observe fast rupture, slow rupture, and the transitions between these regimes.
A. The influence of front type on the loading curve The event shown in Fig. 3 is a fast-slow-fast front. The same simulation also produces fronts that are fast across the entire interface. The associated drops in the loading curve F T /F N have approximately the same amplitude (see Fig. 2). This indicates that the details of the front propagation do not influence the loading curve strongly, at least not for these parameters. Nevertheless, the slow front propagation does have a signature in the loading curve that appears when we zoom in on a few events as in Fig. 2. Namely, the drop in the loading force has a significant change in slope while the slow front lasts, which distinguishes it visually from the force drop associated with a fast-only event.
Let us consider the evolution of F T in more detail. As long as the slider remains pinned, F T increases with the motion of the driving stage (the driving stage moves the end of the driving spring that is not attached to the slider). F T decreases only when the point on the slider where the driving spring attaches, the trailing edge, moves away from the driving stage. This occurs in two distinct ways. First, the trailing edge moves away from the driving stage when the slider deforms in compression during the passage of a rupture front, which happens both for partial slip and full sliding events. Second, the trailing edge moves with the rest of the slider when the entire interface is slipping in a full sliding event. For full sliding events we define the boundary between rupture front passage and full sliding as the moment the rupture front reaches the leading edge, seen e.g. in Fig. 3. The relative amplitude of the F T reductions associated with each of these two motions of the trailing edge depends on the relative stiffnesses of the slider and the driving spring. As seen by the relative amplitudes of the force drops in Fig. 2, with the present parameters the (trailing edge) slip associated with the slider deformation accounts for about a fifth of the net slip of full sliding events.
Further quantification of this feature is presented in Fig. 5, where the force drops associated with either the entire event or with the rupture front passage only are presented for all events during the developed stick-slip regime. We find that the loading force drop occurring during the rupture front passage, regardless of the type of event, is always about four times smaller than the one associated with sliding. Because block motion during full sliding events accounts for the largest part of the net block motion, fast-slow-fast and fast-only events cannot be distinguished from the amplitude of their net force drop. type. We include partial slip events and full sliding events occurring between t = 0.63 s and t = 1.50 s (the developed stick-slip regime). The drops occurring during the passage of the rupture front have comparable amplitude between the partial slip events, fast-slow-fast events and fast-only events. Also, the net force drop in full sliding events have comparable amplitudes for fast-slow-fast and fast-only events.

B. Fast slip and fast front speeds are inertial
In this section we demonstrate that both the fast slip part of the block slip evolution and the fast front speed are of inertial origin. That is, like the bulk wave speeds, these speeds scale as ρ −1/2 , where ρ = M/(LBH) is the mass density of the system.
To isolate the effect of inertia from the stress and frictional state at the interface we have performed four simulations starting from the same state (Appendix A), but with ρ decreased to 1, 1/2, 1/4 and 1/8 of its value in Table I; we changed the mass and kept the system size constant. Figure 6a shows the fast slip dynamics for four neighboring blocks located within the part of the interface where the front speed was high as the front passed, for each of the four simulations. The figure clearly shows that the fast slip speed was modified by the change of density. Figure 6b demonstrates that a ρ −1/2 scaling collapses the data. The reference time t rup for each block is the time of block rupture as defined in Fig. 3, and the block slip is measured with respect to the block position at t rup . Figure 6c shows the front speed as a function of position along the interface for the four simulations, while    6d shows the same data with the front speed rescaled by ρ −1/2 . The fast front speeds are collapsed onto each other by this rescaling. The slow front speeds, in contrast, remain the same in all four simulations (Fig. 6c), and are split from each other by the scaling. This indicates a non-inertial origin of the slow front speed.

C. Predictive power of the front type map
It will be useful to define the distribution φ(f T ) of the forces f T in the junctions attached to one block, and to consider the relationship between φ(f T ) and the effective static friction coefficient of the block. As shown numerically in e.g. [36] and theoretically in e.g. [54,59,70] and Appendix C, for the rupture of heterogeneous systems in which a number of junctions are loaded in parallel, the maximum load that an interface can bear is related to the width of the load distribution and/or threshold distribution of the various junctions. Homogeneous systems (vanishing distribution width) have the maximum possible macroscopic rupture threshold because all junctions will contribute with their maximum force when collective rupture is reached. In contrast, in heterogeneous systems (finite width) the weaker and/or initially more highly loaded junctions will break first, so that when macroscopic rupture occurs, only a fraction of the initial population of junctions will contribute to the total force.
One therefore expects that increasing (decreasing) the width σ of the junction force distribution to make the interface weaker (stronger) favors fast (slow) front propagation. It also makes intuitive sense that higher prestress, hence smaller distance to the breaking threshold, would favor fast front propagation. We have performed simulations where these two parameters are varied systematically while remaining homogeneous along the interface (see Appendix A 3), and the observed front types are presented in Fig. 7. In the arresting region, the fronts are partial slip events, that is, they stop before reaching the leading edge (some of them stop early, some almost reach the leading edge). The region labeled slow includes all those events that have a slow front part, even if the event is fast along most of the interface. They share the characteristic that the events would arrest in the absence of the slow slip mechanism. In the fast region, events are fast across the entire interface.
In Section V we will require series of simulations that stay within the fast-only or the fast-slow-fast front regimes over a range of variation inτ 0 or σ. The front type map makes it straightforward to choose the initial conditions for these simulations. Our choices are indicated with circular and square markers in Fig. 7.
We emphasize that while the arresting, slow and fast regions of the rupture fronts are robust in their relative positions (the fast front region is found at higher values of normalized prestressτ 0 and junction distribution width σ than the slow front region, which is itself found at higher values than the arrest region), the precise locations of the boundaries between them depend also on how the events are triggered (all the simulations in Fig. 7 have the same settings in the triggering region). We also wish to note that the front type map is not a local measure. For example, as mentioned above, many of the events that populate the slow front region are of the fast-slowfast type, meaning that even though the prestress and junction force distributions are the same for all the blocks in the propagation region, the rupture passes some of them as a slow front and others as a fast front.
If these limitations are kept in mind, however, the map is a powerful tool for guiding our intuition. For example, in an interface where there is a sudden change of initial conditions within the propagation region, the map tells us what change to expect in the rupture front. In Initial states in the forbidden region would have had some junctions stretched beyond their breaking threshold and are therefore excluded. The symbols are for simulations used in the following figures: Fig. 8: black star and arrows. Fig. 10: open circles (vertical, only three of the points are shown here). Fig. 11: black filled circles, magenta (gray) filled squares, and multicolored (gray) filled circles (horizontal). , this can lead to a transition from slow to fast front propagation right at the location where the change occurs. Conversely, some (but not all) simulations with a change of opposite sign (from a high to a low prestress, for example), lead to a fast-slow transition.

V. FRONT SPEED RESULTS
In this section we consider the difference in speed between two events that are either both fast, or both slow. In the language of fracture mechanics, rupture occurs when the energy available at the front tip reaches the energy required to break the contacts there. Thus, the speed at which rupture propagates depends on the en- ergy that is already present (related to the stress state), the energy level needed to break the contacts (related to the local strength) and on how quickly the missing energy can be supplied (which depends on the slip motion behind the tip and is therefore transient). Although it has been shown that fracture mechanics satisfactorily describes fast shear rupture events [14,23], in the model, it is more straightforward to argue in terms of forces rather than energies; however, we will see that considerations similar to those in fracture mechanics apply.

A. Front speed is transient
If the front speed did not have transients, or if the transients were short compared to the length scale at which we study the front propagation, the rupture front speed would be uniquely determined by the local state of the interface at the rupture tip, i.e. local stresses, strength, stiffnesses etc. Indeed, various local interfacial parameters have been shown to be correlated with the local front speed (see e.g. [9,17,34]). However, our simulations show that the speed at any point depends not only on the state at that point, but also on the region the front has just passed through. Near the trailing edge, within the length required for the speed of the newly nucleated front to converge, the region behind the front that influences the propagation extends back to the region where the front was triggered, and so in addition to the state of the interface, the precise way rupture starts also is part of what determines front speed.
To illustrate the transient nature of the front speed, Fig. 9a shows the front speed as a function of position for two simulations in which the local stress state and local frictional strength are homogeneous along the propagation region (see details about how these single event simulations were set up in Appendix A 3). We first consider the shorter of the two systems (magenta line with dotted markers). We observe that the front speed changes throughout, that is, even though the local state is the same along the slider, the front speed is not. Note that this is a fast-only front; even though the front speed starts low, its propagation does not depend on the slow front mechanism.
To investigate the convergence length of the transients and the spatial extent over which edge effects dominate we performed another simulation. This simulation differs only in the total length of the sample, which we increased by a factor of five. This gives the black line with circular markers in Fig. 9a. We observe that the two rupture events have very similar front speeds. When the front tip approaches the leading edge, there is a change in curvature that we interpret as an edge effect. This interpretation is supported by the fact that the effect is present near the leading edge of the sample for both the short and the long system, but the transient effect observed at the end of the short sample, around blocks 50-57, is not observed at blocks 50-57 in the longer sample, where these blocks are far from the edge. Figure 9b, where front speed is plotted against inverse position along the interface, shows that the front speed does converge to a finite value. Note that we do not expect a 1/x behavior to hold in general; rather, we plot front speed versus inverse position because it is the simplest function that brings very large x values close to the ordinate axis, and thus facilitates a crude extrapolation to the expected front speed reached in very long samples. Figure 9c and d also indicate convergence. They show two events. The local stresses and strengths are the same in both, but the initialization of the events is different. line with dotted markers: front triggered by driving from the trailing edge, which modifies the stress in the loading region prior to rupture. In all panels, rupture velocity is defined as in Fig. 3. In (b) and (d), the speed for the ten blocks closest to the leading edge are excluded from the plot, to avoid the region where the edge effects dominate.
This leads to front speeds that are initially also different, but which converge to the same value.

B. Front speed depends on local stress state
While we and others [34] argue that the ratio τ 0 /p of shear to normal stress at the interface as the rupture front begins provides insufficient information for predicting the front speed on its own, it remains true within the model, as is also seen in experiments [17], that the front speed has a strong dependence on τ 0 /p. Higher prestress τ 0 results in higher front propagation speed. This result agrees with our intuition: for a given strength, a region that is more highly prestressed requires less stress change to start slipping, and it releases more energy when doing so. Both favor higher front speed. This is also consistent with Fig. 7, which showed that increasing the initial shear stress in the propagation region will bring a system into a regime where fronts are fast across the entire interface. Figure 10a shows the front speed for events in which the initial stress state was varied systematically while remaining homogeneous along the interface (a similar variation in the junction force distributions' width σ is the topic of Section V C). In agreement with our argument above we observe in Fig. 10a that higher prestress (dark red line) corresponds to higher mean front velocities.
Note from Fig. 10 that in the model the longitudinal wave speed c L = 1677 m/s does not set an upper limit for the front propagation speed. Other studies of bi-material interfaces, both experimental (see e.g. [71]) and numerical (see e.g. [34]), also found that when the substrate is stiffer than the slider, as is the case here, the rupture front speed can exceed the longitudinal wave speed of the slider.

11
C. Front speed depends on local strength In the same way as for the influence of stress state on front speed discussed above, it makes intuitive sense that for a given prestress, the front propagates faster if the frictional strength, i.e. the threshold stress for slip inception, is lower. In the model, the frictional strength is controlled to first order by the values of f thres and f slip on the junction level, but even when these are the same everywhere, the local strength can vary. We show in Fig. 18 (Appendix C) how increasing the width, σ, of the initial junction force distribution, φ(f T ), results in a lower effective static friction threshold, and we expect that if conditions are otherwise the same, this will result in faster front propagation. In the case (not studied here) of a more complete model, variations in local strength would also be due to variations in individual strength of the junctions attached to each block.
We observe in Fig. 11a that for the same prestress, higher σ does result in higher front speed. Within our sample size, which is the same as in previous sections and was selected from experimental parameters, the front speed is changing throughout the front propagation. In order to investigate the effect of σ quantitatively and to isolate it from the effect of prestress, we compare each event to a reference event having the same prestress state and with width σ = 0. We measure along the entire interface and find average and min/max values. The results are presented in Fig. 11b. By rescaling the front speeds with the reference speed for σ = 0 the similarity between the two series for different prestress becomes apparent. The range of results for each value of σ due to the transient nature of the front speed is large, but there is a clear trend that higher σ (weaker interfaces) correspond to higher front speeds. The effect of σ on speed can be very significant, giving up to a 50% relative increase for the data shown in Fig. 11b. Figure 11a and b show data for fast-only fronts. The influence of σ on front speed is qualitatively similar for slow fronts. This is shown in Fig. 11c and d. For the fronts which have both fast and slow propagation parts, we used visual inspection of the underlying data in Appendix D to determine which part of the fronts to include. Our selection can be read off of the horizontal extent of the lines in Fig. 11c. Again, the effect of σ on the front speed is significant, with an increase of about one order of magnitude when σ/(f thres − f slip ) is increased from 0.12 to 0.47.
Note that the ratio σ/(f thres − f slip ) in interfaces that support front propagation remains small even if f thres − f slip is chosen to be small or zero. This is because σ, the width of the distribution forces in a block's junctions, also vanishes in this limit. The explanation is twofold. First, the average junction force prior to rupture must be equal to or larger than f slip , otherwise the front arrests. Second, the force in each junction cannot exceed f thres . Thus the width of the junction force distribution vanishes as f thres approaches f slip . The motion of the rupture front tip and the motion of the material of the slider are interrelated: the slider cannot move while the interface is in the pinned state (the slip depends on the front), and the rupture propagates as deformation of the slider transfers stress and energy to the front tip (the front depends on the slip). In [36] we demonstrated that the slow front speed in the model is proportional to the slow slip speed and worked out the constant of proportionality. In this section we show that the fast front speed in the model depends on the fast slip speed through the exact same relationship. Figure 12a shows the fast front speed vs fast slip speed for the blocks between x = 3 cm and x = 13 cm in a series of simulations with controlled initial states and variation in mass density, normal force, junction distribution width, prestress and bulk and interface stiffnesses (Fig. 12c shows the same data on logarithmic axes). These blocks are all in the propagation region and away from the triggering zone and the leading edge. The front speed at each point along the interface is measured in the same way as in Fig. 3. The corresponding slip speed is the average block slip speed measured from the time block i starts slipping to the time when the block 5 mm away (block i + 2) starts slipping. We choose this time interval because it is during this period that the motion of block i most directly affects the rupture front propagation. Figure 12b shows the data in Fig. 12a with the fast slip speed in each simulation rescaled according to the parameters used (Fig. 12d shows the same data on logarithmic axes). We also include (in black) the slow front and slow slip data from [36, Fig. 3D]. Arguments for the scaling were presented in [36, SI Equations]; we repeat only the conclusions here. The rescaled slip speed is v slip, rescaled = v slip kil0 τ thres −τ0 , where k i = j k ij is the stiffness of the connection between a block and the interface, l 0 = 7 mm is the characteristic decay length of the shear stress field, which depends on the bulk to interfacial stiffness ratio k/k i , τ thres = N s f thres , and τ 0 is the initial shear stress before the event. The normal force enters in the scaling indirectly because it modifies f thres . We expect deviations from the rescaling due to two sources that have not been included in the scaling equation. First, the effective force threshold τ eff thres on a block with a given junction force distribution is always less than τ thres , as discussed in Appendix C. The data where the junction force distribution was varied collapses better when τ eff thres replaces τ thres in the scaling relation. However, because τ eff thres depends on the microscopic state of the interface, we choose to omit this correction in order to keep the factors in the scaling equation at the mesoscopic level. Second, the shape of the shear stress field also appears in the argument for the scaling, but not in our final scaling equation, which only includes the characteristic stress decay length l 0 ; we have kept its variation between simulations to a minimum by keeping the ratio  Fig. 7 as an arbitrary reference, we have varied the junction force distribution φ(f T ) (magenta/gray crosses), the prestressτ 0 (green/gray dots), the mass density ρ (blue / dark gray squares), the bulk and interfacial stiffnesses (keeping k/k i constant) and ρ (cyan / light gray circles), and the normal force F N and ρ (red/gray diamonds). k/k i constant. Figure 12b shows that the above rescaling allows all data from Fig. 12a to nicely collapse on a single straight line, going through the origin and having a slope not far from 1. This demonstrates that the front speed is directly proportional to the concurrent slip speed, with the proportionality coefficient being kil0 τ thres −τ0 . Strikingly, the data for slow slip and slow fronts collapse on the very same line, as clearly seen in Fig. 12d. Figure 13 shows the rupture front speed vs prestress ratio τ /p for several events where the shear prestress and the junction force distributions were varied systematically. In order to both enable the tuning of prestress and force distributions and to mimic the experimental driv-13 ing conditions more closely, as in the full simulations, initial conditions were prepared in two steps. First, homogeneous distributions and prestresses were assigned. Then, the trailing edge driving spring was allowed to move, eventually triggering the front. This raises the prestress most prominently near the trailing edge, but for the present parameters (which are found from comparison with experiments) the ratio of shear stress profile decay length, l 0 , to system length is such that the shear stress is non-negligibly changed everywhere. For each event, rupture speed is measured at all points between x = 2.5 cm and x = 11 cm, i.e. avoiding the immediate vicinity of the system edges. Figure 13 is analogous to Fig. 3 in Ben-David et al. [17], which shows the rupture front speed in selected points away from the sample edge vs prestress. Both figures show the same trend, namely that the front speed increases with increasing prestress, and that there exists a continuum of front speeds ranging from close to zero (slow) up to super-shear wave speeds. In addition, the numerical data indicates a systematic effect of varying the width of the junction force distribution, in a way which is completely consistent with the observations reported in section V C: larger widths yield lower local strength and thus promote faster fronts.

VI. HISTORY DEPENDENCE OF MESOSCOPIC STATIC FRICTION COEFFICIENT
In the previous section we argued that the rupture front speed depends on the local strength, i.e. the effective static friction threshold on the block level.
Here we investigate the extent to which variations in the mesoscale strength occur even when the individual micro-junctions attached to the block all have the same strength. This effect comes in addition to spatial heterogeneities in stresses, due for instance to previous ruptures of the interface, which also cause spatial variations in the mesoscale strength.
We showed in [59] (and recall in Appendix C) how the local strength depends on the distribution of forces among the micro-junctions of a block. We also showed how, for simple velocity profiles, the distribution of forces at block arrest depends on the deceleration that led to this arrest. More fundamentally (in the model), the distribution of forces at arrest depends on the slip profile x(t) of the block during the re-pinning of its junctions. As simplified velocity (and thus slip) profiles were studied in [59], we focus here on the more complicated spontaneously occurring events in the full simulations. We are primarily interested in the effective static friction threshold µ eff s = τ max /p after an event, where τ max is the largest tangential force that the block will bear before it starts to slip in the next event. µ eff s can be calculated directly from the junction force distribution if we assume that the rupture is fast compared to the mean re-pinning time. As we will show, however, µ eff s can also be estimated from the preceding slip dynamics.

A. Gini coefficient: definition and measurement
The junction force distribution φ(f T ) of a block evolves with the slip history of the block [59]. As the block starts moving, slips and comes to rest, junctions are continually breaking and reforming, and the complete distribution φ(f T ) depends on the full slip velocity history v(t). In practice, however, it is the very end of the preceding Appendix B) how ibution of stretchblock. We also les, the distributhe decelerations. he distribution of profile x(t) of the ctions. As simplire studied in [20], ccurring events in rested in the e↵eccan be calculated assuming that the repinning time.
Gini coe cient as plicated block slip at the Gini coed predictor of the block will have in d measurement law that we studtribution (f T ) of the block. As the to rest, junctions , and the complete l slip velocity hishe very end of the fluence on , and strength of after block slip motion. etween full sliding strength of . The time. Second, we y that needs to be at the chosen time. kwards from t 0 , we was one junction it is at t 0 ; x(t 0 ) can be neglected oints that existed during t 2 [t 0 , t 0 ]. R , so that [t 0 , t 1 ] occurring after t 0 .
blocks remain at int can be chosen. ases where [t 0 , t 1 ] events that set up simulations this is ip history. Third, compare it to slip lly or with simple its influence on otion x(t). This We showed in [20] (and repeat in Appendix B) how he local strength depends on the distribution of stretchngs among the microjunctions of a block. We also howed how, for simple velocity profiles, the distribuion of stretchings at arrest depends on the decelerations. ore fundamentally (in the model), the distribution of tretchings at arrest depends on the slip profile x(t) of the lock during the repinning of the junctions. As simplied velocity (and thus slip) profiles were studied in [20], e focus here on the spontaneously occurring events in ull simulations. We are primarily interested in the e↵ecive static friction thresholdμ e↵ s , which can be calculated rom the stretching distribution when assuming that the upture is fast compared to the mean repinning time.
In Section VI A we introduce the Gini coe cient as scalar, integrated quantifier of complicated block slip ynamics. In Section VI B we show that the Gini coeient of a block's slip motion is a good predictor of the ↵ective static friction threshold the block will have in ts next event.

A. Gini coe cient: definition and measurement
A feature of our microscopic junction law that we studed in [20] is how the junction force distribution (f T ) of block evolves with the slip history of the block. As the lock starts moving, slides and comes to rest, junctions re continually breaking and reforming, and the complete istribution (f T ) depends on the full slip velocity hisory v(t). In practice, however, it is the very end of the receding event that has the largest influence on , and his allows us to predict the e↵ective strength of after n event from a characterization of the block slip motion. First, we pick a point t 0 in time (between full sliding vents) at which we wish to predict the strength of . The locks should be at rest at this point in time. Second, we dentify the segment of the slip history that needs to be onsidered to predict the strength of at the chosen time. ollowing each block's motion x(t) backwards from t 0 , we efine t 0 as the time when the block was one junction reaking length s m away from where it is at t 0 ; x(t 0 ) (t 0 ) = s m . The slip history before t 0 can be neglected ecause all the junction attachment points that existed efore t 0 have been broken and renewed during t 2 [t 0 , t 0 ]. urther, we define t 1 = t 0 +t R + t R , so that [t 0 , t 1 ] ncludes most of the junction renewal occurring after t 0 . e require t 1 < t 0 , but as long as the blocks remain at est, any t 1 that satisfies this constraint can be chosen. t is possible to create pathological cases where [t 0 , t 1 ] ails to include the bulk of repinning events that set up (f T ) at t 0 , but in practice, in all our simulations this is he most important segment of the slip history. Third, e characterize v(t) for t 2 [t 0 , t 1 ] and compare it to slip rofiles that can be treated analytically or with simple umerical integration. The feature of v(t) that dominates its influence on s how evenly it distributes the slip motion x(t). This makes intuitive sense: If the block comes to rest or nearly to rest during a time interval much shorter thant R , most of the junctions will be renewed while the block is at rest, within a small interval of stretching, and (f T ) will be narrow. If the block moves uniformly, the junctions will distribute more uniformly. We assume here that the rate of junction renewal is constant, a simplification that is a reasonable approximation when the block comes to rest after full sliding. A robust measure of inequality is the Gini coe cient, which we introduce in Fig. 16. It is commonly used to characterize the inequality of the income or wealth distribution in a population [43], but can be applied to our case without modification. The Gini coe cient of a set is defined from the areas in FIG. 16 as G = A/(A + B) (which equals 2A, since A + B = 1/2). B is the area under the cumulative of the share assigned to each element in the set when the elements are ordered by the size of their share. A is the area between the cumulative and line of equality. The line of equality is the cumulative for a set where each element has the same share. Thus, if each element has the same share, G = 0, while in the other extreme where one element has everything and the other elements have zero share, G = 1. To measure the Gini coe cient of the block slip during t 2 [t 0 , t 1 ] we create a set consisting of N e = 10 elements. We define t elem = (t 1 t 0 )/N e and assign x share to each element so that the elements tile the time interval. That is, the share of the first element is the block slip that occurred between t 0 and t 0 + t elem , the share of the second element is the block slip that occurred during [t 0 + t elem , t 0 +2 t elem ], and so on. Equivalently, each element's share is the average slip velocity during the element's time interval (net makes intuitive sense: If the block comes to rest or nearly to rest during a time interval much shorter thant R , most of the junctions will be renewed while the block is at rest, within a small interval of stretching, and (f T ) will be narrow. If the block moves uniformly, the junctions will distribute more uniformly. We assume here that the rate of junction renewal is constant, a simplification that is a reasonable approximation when the block comes to rest after full sliding. A robust measure of inequality is the Gini coe cient, which we introduce in Fig. 16. It is commonly used to characterize the inequality of the income or wealth distribution in a population [43], but can be applied to our case without modification. The Gini coe cient of a set is defined from the areas in FIG. 16 as G = A/(A + B) (which equals 2A, since A + B = 1/2). B is the area under the cumulative of the share assigned to each element in the set when the elements are ordered by the size of their share. A is the area between the cumulative and line of equality. The line of equality is the cumulative for a set where each element has the same share. Thus, if each element has the same share, G = 0, while in the other extreme where one element has everything and the other elements have zero share, G = 1. To measure the Gini coe cient of the block slip during t 2 [t 0 , t 1 ] we create a set consisting of N e = 10 elements. We define t elem = (t 1 t 0 )/N e and assign x share to each element so that the elements tile the time interval. That is, the share of the first element is the block slip that occurred between t 0 and t 0 + t elem , the share of the second element is the block slip that occurred during [t 0 + t elem , t 0 +2 t elem ], and so on. Equivalently, each element's share is the average slip velocity during the element's time interval (net FIG. 14: The areas that define the Gini coefficient, a measure of inequality. event that has the largest influence on φ, and this allows us to predict the strength of φ after an event (which is the effective static friction threshold for the next event) from a characterization of the block slip motion. First, we pick a time t between two full sliding events at which we wish to predict the strength of φ. The blocks should be at rest at this point in time. Second, we identify the segment [t 0 , t 1 ] of the slip history that needs to be considered to predict the strength of φ at the chosen time. Following each block's motion x(t) backwards from t , we define t 0 as the time when the block was one junction breaking length s m = f thres /k ij away from where it is at t ; x(t ) − x(t 0 ) = s m . The slip history before t 0 can be neglected because all the junction attachment points that existed before t 0 have been broken and renewed during t ∈ [t 0 , t ]. Further, we define t 1 = t 0 +t R + δt R , so that [t 0 , t 1 ] includes most of the junction renewal occurring after t 0 . We require t 1 < t , but as long as the blocks remain at rest, any t 1 that satisfies this constraint can be chosen. It is possible to create pathological cases where [t 0 , t 1 ] fails to include the bulk of re-pinning events that set up φ(f T ) at t , but in practice, in all our simulations this is the most important segment of the slip history. Third, we characterize v(t) for t ∈ [t 0 , t 1 ] and compare it to velocity profiles that can be treated analytically or with simple numerical integration.
The feature of v(t) that has the most prominent impact on φ is how evenly the slip motion x(t) is distributed in time. If the block comes fully or nearly fully to rest during a time interval much shorter thant R , most of the junctions will be renewed while the block is at rest; they will therefore have approximately the same stretching and φ(f T ) will be narrow. If the block moves more uniformly as it comes to rest, the stretching of the junctions will distribute more uniformly. This simple argument assumes that the rate of junction renewal is constant, a simplification that is a reasonable approximation when the block comes to rest after significant slip motion [59].
A robust measure of inequality is the Gini coefficient, which we introduce in Fig. 14. It is commonly used to characterize the inequality of the income or wealth dis-tribution in a population [72,73], but can be applied to our case without modification. The Gini coefficient of a set is defined from the areas in Fig. 14 as G = A/(A + B) (which equals 2A, since A + B = 1/2). For any set of elements, each having a share of a some quantity (e.g. wealth), the line separating A and B is constructed as follows. Elements are first ordered by the size of their share and placed along the abscissa. The ordinate for a given share then corresponds to the sum of all shares that are smaller than the chosen share. The line thus defines the cumulative distribution of shares. B is the area under the cumulative and A is the area between the cumulative and the line of equality. The line of equality is the cumulative for a set where each element has the same share. Thus, if each element has the same share, G = 0, while in the other extreme where one element has everything and the other elements have zero share, G = 1. To measure the Gini coefficient of the block slip during t ∈ [t 0 , t 1 ] we create a set consisting of N e = 10 elements. We define ∆t elem = (t 1 − t 0 )/N e and assign x share to each element so that the elements tile the time interval. That is, the share of the first element is the block slip that occurred between t 0 and t 0 + ∆t elem , the share of the second element is the block slip that occurred during [t 0 + ∆t elem , t 0 + 2∆t elem ], and so on. Equivalently, each element's share is the average slip velocity during the element's time interval (net slip equals average slip velocity times the length of the time interval, but the common normalization factor ∆t elem can be ignored). Figure 15a shows block slip profiles x(t) for constant deceleration slip. These slip profiles have the advantage that they are completely characterized by the deceleration amplitude; they help build intuition before we consider the full simulation data, where none of the basic kinematic quantities (slip, velocity, acceleration) are even approximately constant. In Fig. 15a the time interval has been shifted so that t 0 = 0. With the parameters in Table I, (t 1 − t 0 )/t R = 1.3. Figure 15b shows the effective static friction, i.e. the strength of φ(f T ), resulting from the slip profiles in Fig. 15a. We showed in [59, Fig. 9] that for (i) initial velocities large enough that the blocks move longer than the junction breaking length s m before coming to rest, and (ii) reasonable initial junction distributions, the effective static friction for constant deceleration slip depends only on the deceleration value and the junction law parameters; here we use the same initial conditions as in [59, Fig. 9] and the junction law parameters in Table I. Figure 15c shows the Gini coefficient of the slip profiles in Fig. 15a as a function of their constant deceleration parameter (Gini coefficients were calculated using [74]). The shape of the curve can be intuitively understood. For deceleration close to zero, the slip motion from t 0 to t 1 occurs with low and nearly constant speed, so the cumulative is close to the line of equality, and G ≈ 0. For large deceleration values, the block comes to rest sharply, [t 0 , t 1 ] includes a long period where the block is at rest, and the cumulative is more peaked. In the limit of very large deceleration values, G → 1. We can combine Fig. 15b and Fig. 15c to obtain a prediction for µ eff s vs G. This gives the full drawn lines in Fig. 17. Figure 16 shows example slip profiles from full sliding events in full simulations. Two events are shown, with the data in c being a detailed view of the data in a and the data in d being a detailed view of the data in b.

B. Gini coefficient predicts effective static friction
For the limiting cases of a block coming to rest very abruptly or very gradually, the Gini coefficient and the effective static friction threshold can both be predicted exactly from simple arguments. The abrupt stopping case is the limit where the block comes to rest from a speed high compared to s m /t R within a distance short compared to s m and a time short compared tot R . In this case the entire motion between t 0 and t 1 will be assigned to a single element and the Gini coefficient will be G = 1. The distribution of junction forces will be a δ-function as all junctions reform after the block has stopped moving, and so τ max /p = f thres /f N and Colored line with markers: a moving average filter of width 1 ms (=t R /2) was applied to x(t) and then N e = 10 equally sized intervals were defined (markers indicate intervals' boundaries). The smoothing helps avoid negative increments that while possible to include are not part of the basic formulation of the Gini coefficient.
the rescaled effective static friction threshold isμ eff The gradual stopping case is the limit where the block spends a long time at velocities small compared to s m /t R . In this case the velocity is nearly constant from t 0 to t 1 and the Gini coefficient will be G = 0. A uniform junction force distribution is set up (see [59, Section III A] for a detailed account). From equation (C2) we find µ eff s = (f thres − f slip )/(2f thres ). As a slipping block slows down the slow slip mechanism becomes important. If the amount of slow slip is large compared to the junction breaking length s m , the motion is always far from the abrupt stopping limit. Consequently, to span out the range of Gini coefficients from 0 to 1 we have performed simulations where we in- directly varied the relative amplitude of slow slip and s m by varying f thres and f slip . Figure 17 shows results from this series of full simulations. As expected, the variation between events within one simulation is small compared to the variation between simulations. In all cases, the Gini coefficient is a good predictor of effective static friction. In particular, we find near equality between G and µ eff s for G > 0.5. The agreement between the line corresponding to constant deceleration slip and the simulation data is rather good, but far from perfect, as expected given the approximations inherent in describing the complicated block slip motion arising in the full simulations with a single number. Important differences in the underlying motion are: (i) In the full simulations the duration of sliding events is on the order oft R , see Fig. 16. To reach steady state junction distributions, as was assumed in [59], a longer sliding period is necessary. (ii) In the full simulations, the oscillations in the block velocity set up junction force distributions that are non-smooth. This is why the simulation data is not bounded by the same lower limit as the constant deceleration data.

VII. DISCUSSION
We begin this section with an overview of the robustness of the observed front dynamics to variations in model parameters. We then interpret our results and compare them to earlier work in the literature.

A. Robustness
To elucidate the robustness of the system dynamics to variations in the model parameters we have performed simulations where one or two parameters at a time were varied from their values in Table I. We find that most of the parameters are in a range where the results remain qualitatively the same under halving and doubling of the parameter value. Others, notably the bulk and interface compliances, are in a range where these changes lead to significant changes in the system dynamics.
As mentioned in Section II, T (t R ) is a simplified way of modeling the distribution of times after which slipping micro-junctions relax. It is therefore reassuring that the front dynamics are not sensitive to the details of T (t R ): (i) We observe fast-slow-fast fronts in a simulation with an exponential T (t R ). (ii) We varied the average, t R , systematically and observed that this changes the speed, but not the nature of slow slip and slow fronts. (iii) We have performed full simulations where we halved and doubled the width, δt R , from its value in Table I. Although the details of the fronts change, we still observe precursors, fast-only events, fast-slow-fast events and intermediate partial slip events.
The number of junctions per block, N s , could be adjusted to match the area density of junctions in a given experiment. However, for N s 1 we expect to recover the results from [59], where instead of tracking individual junctions we described the state evolution in terms of continuum distributions of pinned and slipping junctions. With N s = 100, the present simulations appear to be within this continuum limit: when we compare to full simulations with N s = 50, 200 and 500 we again observe precursors, fast-only events, fast-slow-fast events and intermediate partial slip events, and although due to the randomness in the junction dynamics the events are not exactly the same between these simulations, they are very similar.
There are four stiffness parameters in the model: the driving spring modulus K; the elastic foundation modulus k f ; Young's modulus E, from which the bulk spring stiffness k follows; and the interface spring stiffness k ij . We have performed simulations where the driving spring modulus K was increased twofold, tenfold and hundredfold from its value in Table I. This changes the relatively regular stick-slip of in the loading curve to a more chaotic stick-slip pattern, and as predicted in Section IV A it causes the amplitude of the drops in F T during partial slip events to increase. We still observe both fast-only and fast-slow-fast fronts. For the elastic foundation modulus, k f , doubling and halving its value in full simulations still allows the full range of system dynamics. For k f k, the normal stress profile p(x) is essentially flat, as significant variations in p would then require large vertical deformations of the slider, and the experimental stress profiles of [6,17] will no longer be reproduced.
For the bulk and interface stiffnesses k and k ij we have shown in Fig. 12 that when their ratio remains constant, changes to the value of these parameters affect the slip and front speed, but dynamics remain qualitatively the same. When their ratio changes, however, qualitative changes to the system dynamics occur. Doubling k or halving k ij increases the loading region where the interface stress is significantly influenced by the increasing force in the driving spring between events. This leads to multiple front tips at different locations propagating simultaneously instead of one front propagating from the trailing to the leading edge, so that a single front speed cannot readily be defined. Halving k or doubling k ij , on the other hand, reduces the extent of the loading region, so fronts can still be defined as before, but these settings favor fast front propagation, so very few fast-slow-fast fronts were observed. Of course, the slow slip mechanism is still active, but it is masked by the fast slip and fast front propagation.
Variations in the force levels in the junction evolution law, f thres and f slip , are included in the scaling results explained in Appendix C. We also expect that parameter values can be found for which full simulations produce precursors, fast-only events, fast-slow-fast events and intermediate partial slip events with larger values of f thres and f slip than we use. Because this change modifies both the triggering and propagation of events, other parameters would have to be changed simultaneously.

B. Interpretations and comparisons
In Section IV A we identified a signature of slow front propagation in the macroscopic loading curve. If this signature turns out to also be present in the experiments where slow fronts are known to occur, it could be a useful first indicator of slow fronts in experiments where the nature of rupture front propagation is unknown. If the reduction in F T during front propagation is sufficiently large to be reliably measured, a fast-slow transition may show up in F T (t) as a reduction in slope. If the front transitions back to fast propagation or reaches the leading edge, a change back to large negative slope may occur in F T (t).
In Section IV B we demonstrated by changing the slider mass density ρ that, like the bulk wave speeds, the fast slip and fast fronts in the model are of inertial origin. As, in the model, fast slip and slow slip are due to distinct mechanisms, it follows that the ratio of their speeds can vary if system parameters are changed. For instance, in the model: fast slip speed is independent of and slow slip speed is inversely proportional to the time scalet R which determines interface healing after rupture; fast slip speed is inversely proportional the square root of slider mass density and slow slip speed is independent of slider mass density. In contrast, reducing (increasing) the interface stiffness increases (reduces) both slow and fast slip speeds without affecting their ratio. Because front speed is proportional to slip speed across the velocity scales that we have observed, see Fig. 12, it follows that large variations in the ratio of slow to fast front speeds is also possible between systems. This paradoxically suggests that the usual definition of slow fronts in terms of their velocity (one to several orders of magnitude slower than the Rayleigh wave speed), although natural in the system in which they were discovered [3], may not be satisfying. In our model, and possibly in other systems, a definition based only on fronts' speed does not capture the physical characteristics of slow fronts -dynamic fronts that in the model propagate due to the intrinsic relaxation dynamics of the frictional interface after arrestor differentiate them from other front types. This definition may lead to misinterpretations in systems where the velocity range of slow fronts significantly overlap the velocity range of fast fronts (of the order of sound speed) or quasi-static fronts (speed proportional to the rate of external loading).
In Section IV C we showed that the type of front observed (arresting, fast-slow-fast, or fast-only) depends systematically on both the strength and the prestress of the interface, with weaker interfaces and higher prestresses favoring faster front propagation. We then illustrated that although the front type map in Fig. 7 is not a local measure, because front speed is transient and also depends on the interface state behind the front tip, the map can still give useful predictions of the way front type changes with sudden changes in interface conditions. This could help in the interpretation of experiments on interfaces with heterogeneous friction properties like [75], where material heterogeneities in the form of rock pebbles embedded at a gel-glass interface were observed to cause the rupture front speed to sometimes increase, sometimes decrease.
In Section V B we showed that, in agreement with experiments and other numerical work, front speed increases with increasing shear to normal prestress ratio τ /p. However, knowing the local τ /p alone is insufficient information to predict the front speed. Indeed, in Section V A we showed that in the model and for experimentally realistic system sizes, front speed is transient at all points and as such the selection of its value is both history-dependent and non-local. In Section V C we showed that both fast and slow front speeds also depend on the local strength of the interface, which is generally out of experimental reach, especially when the interfacial junctions bear some randomness, e.g. in their size or stiffness.
In [36] we demonstrated that the slow front speed in the model is proportional to the slow slip speed and worked out the constant of proportionality. In Section V D we found that the same constant of proportionality relates fast slip to fast front speed. The proportionality follows from relating how far the blocks behind the front tip have to move to advance the rupture, to how quickly they move (the slip speed). Bar Sinai et al. [28, equation (9)] (see also [27, equation (9)]) found a similar result in an analytical model with a rate-and-state-type friction law. Taken together, these model results suggest that the origin of the front speed, be it slow or fast, lies in the relevant slip speed during front propagation, the selection of a fast or a slow front regime being the direct consequence of the selection of a fast (inertial) or a slow (driven by the internal relaxation of an arrested interface) slip regime. We therefore suggest that a way to better understand front type selection and the transition between front types is to investigate the physical origin of slip motion and its time evolution during an event.
In Bar Sinai et al. [28], the slow slip speed was selected as the speed for which the steady-state friction law had a minimum. However, the existence of such a minimum is not a necessary condition for slow fronts to appear in a system. As a matter of fact, in our model, and for the parameters used, the steady state friction law is purely decreasing [36]. The apparent contradiction disappears as soon as one remembers that in our model, the slow slip speed arises from a completely different physical mechanism, namely the time distributed relaxation of the interface after its rupture and arrest.
For the scaling of the shear to normal stress ratio we have usedτ 0 = (τ 0 /p − f slip /f N )/(f thres /f N − f slip /f N ) throughout. It would have been consistent with the observed dependence of front speed on both local stress state and local strength to replace f thres /f N with τ max /p = µ eff s to obtainτ eff 0 = (τ 0 /p−f slip /f N )/(τ max /p− f slip /f N ), and use this in Fig. 10b. Nevertheless, we have chosen to stick with a scaling relation that depends only on input model parameters and not on the emerging effective strength of the interface, because quantitative validation of the more advanced scaling is not directly possible due to the transient nature of the fronts. In an application where fronts were known to have converged, it would be interesting to try theτ eff 0 scaling when plotting the fronts' stationary speed against prestress. Similarly, replacing τ thres with τ max in the scaling relation between front speed and slip speed in Fig. 12 would improve the scaling of the data where φ(f T ) was varied, at the cost of introducing emerging properties into the scaling.
Ben-David et al. [17, Fig. 3] reported a relationship between front speed and shear to normal prestress ratio. In Section V E we showed that although our results are quantitatively different (in particular, the range of τ 0 /p values that are stable in the model is about half of that in the experiments), qualitatively we capture the experimental observation well. Further, our results for front speed transients (Section V A) and dependence on local strength (Section V C) are possible explanations for the relatively large range of front speed values that were observed for each value of τ /p in [17]. Namely, we expect variation in local strength to be present in the experiments, probably to a larger degree than in the simulations, as we have explicitly avoided spatial heterogeneities in the friction parameters. It is also possible that the experimental fronts share with our simulations the property that over the sample length studied, the front speed is in the transient regime: the fronts in e.g. Fig. 2 of [17] actually change their speed throughout the interface, but whether this is due purely to changing interface conditions or also to a "true" transient behavior is difficult to assess.
We argued in the model description (Section II) that a time-dependent rule for the reformation of junctions has experimental justification. In the model, a consequence of this microscopic friction law is that the local microscopic state, specifically the distribution of forces in the junctions, depends on the block slip dynamics of the preceding event [59]. The microscopic junction state in turn determines the effective static friction, which affects the front propagation. In Section VI we showed that the Gini coefficient, an integrated estimator of the non-uniformity of the block slip history, predicts the effective local static friction even in full simulations where the block slip dynamics is highly complicated. We were thus able to link two very different aspects of the system behavior. Further, while the distribution of junction forces and the resulting effective friction are difficult to measure experimentally, the local slip history can more readily be measured, at least on the side of the slider [see e.g. 23,60] or in setups which track the motion of patterns or markers at the interface [see e.g. 4,5]. We recognize that there are challenges that need to be overcome to adjust our protocol for measuring the Gini coefficient and apply it to experiments or real systems like seismic faults. However, if the length scale for junction breaking and the time scale for junction reformation could be determined, we believe that the underlying idea, namely to characterize the crucial part of the slip motion with a robust measure of inequality and to use it to predict the (still unknown) strength of the interface, is applicable.
We have studied the onset of sliding friction in a model that couples a microscopic friction law involving the state of a large number of individual junctions to a 2D elastic solver. In the range of possible combinations of friction and bulk models outlined in the introduction this is, we believe, an example of bringing together a moderately complicated friction law with a moderately complicated bulk law. For example, the friction law excludes both aging and a distribution of junction strengths, and the bulk law, while 2-dimensional, excludes among other things plasticity and melting, and is solving for the deformation of the slider only. Nevertheless, the simulations reproduce many of the features of the spatiotemporal dynamics in experiments to which the model was attuned, and within the model we were able to use the more complete information of the system and the fine control of initial conditions available in simulations to better understand how some of the central material properties and system state characteristics modify the dynamics. We believe important avenues for future research to be the inclusion of the bulk deformations of the substrate, which brings with it the complication of defining frictional properties at a non-planar interface, and the parametrization from more fundamental models of the individual junction evolution laws. The simulations we have performed can be grouped in three categories. (i) Full simulations. In these the sample is initially unstressed; then the normal force is applied; then the shear force is applied through the driving spring, and the spring drives the slider through tens of events: precursors, full sliding events, and possibly partial slip events between the full sliding events. (ii) Simulations that restart at a point within a full simulation. In these we make well-defined changes to identify or highlight a particular mechanism, for example the slow slip mechanism. (iii) Smaller simulations where the shear and normal prestresses are controlled and a single event is studied. In these we have done triggering both through (iii)a breaking the junctions in a predefined region, and (iii)b motion of the driving spring as in the full simulations. In this appendix we explain how each type of simulation was set up, and give parameter values in Table I.

Full simulations
In full simulations, the slider is initialized with full normal load F N and no tangential load F T by gradually applying F N without allowing springs to break, a technicality required because the normal forces on the springs, f N ij , start at zero and therefore springs, if allowed to, would break under any stretching. The unique equilibrium is found through damped relaxation of typical duration 10 ms. After relaxation, we check that no spring is stretched beyond its strength and introduce the driving spring starting from zero applied driving force F T . Then F T increases as the driving point moves to the right with speed V . Full simulations were used for Fig. 2, 3, 4, 5, 16 and 17.

Restarted simulations
It can be difficult to analyze events from the full simulations, because the initial conditions for each event arise dynamically from the preceding events. To make direct comparisons between events possible, restarted simulations are run, which begin from an interesting state arising in a full simulation, but continue with modified model parameters or state properties. Restarted simulations were used for Fig. 6.

Single event simulations
Single event simulations, like restarted simulations, were performed to elucidate the qualitative and quantitative importance of individual parameters or interface states by varying them individually, keeping the rest of the simulation setup unchanged. To simplify analysis, these systematic studies were done with different normal forces and different initialization from the full and restarted simulations. The normal force boundary conditions on the top and bottom were exchanged by setting a constant normal force p i = F N /N x on all blocks i at the interface. To maintain stability against global rotation, the top blocks interacted with an elastic ceiling with the same properties as the elastic foundation used in the full simulations.
To obtain an initial state with a prescribed interfacial shear stress profile we turned the interface springs off during the initialization. In their place we added to each bottom block the force corresponding to the shear stress to be prescribed. We also introduced the driving spring, but let V = 0. During relaxation, the sample moved along the x-axis until the force in the driving spring balanced the net force from the interfacial shear stress. To get rid of oscillations more efficiently we added damping forces −α(ẋ i ) on the motion of all blocks. After relaxation, the extra forces and the extra damping were turned off and the interfacial springs were introduced, with their attachment points x ij chosen such that the net force on each block was unchanged and the desired distribution of spring forces, φ(f T ), appeared. We then waited a few timesteps to ensure that the transition from pre-to postrelaxation involved no force discontinuities.
For type (iii)a single event simulations, instead of driving the system with V = 0 until rupture was triggered, we started fronts by simultaneously de-pinning all junctions for all blocks to the left of x trigger . The shear stress in the triggering region has a strong influence on the rupture fronts, and in order to simplify the comparison  [36]. Parameters above the horizontal line were also used, in the same way, in [9]. of results between simulations we used a constant valuē τ trigger = 0.3. The triggering and propagation regions and the initial stress configuration for this type of simulations are shown in the middle row panels of Fig. 8. Type (iii)a simulations were used for Fig. 7, 8, 9, 10, 11, 12 and 19. For type (iii)b single event simulations we matched the experimental event triggering conditions more closely. Again, homogeneous stress and junction force distribution states were prepared as described above. Then, an event was started by moving the driving spring as in the full simulations. This modifies the shear stress state before the event starts, particularly near the loading point. As in our other simulations, the width of the junction force distribution φ affects the effective static friction threshold µ eff s and by that the energy released as the event starts. To isolate the effect of varying parameters in the propagation region, we kept the width of the junction force distribution constant at σ load /(f thres −f slip ) = 0.143 to the left of x load = 22.5 mm. Type (iii)b simulations were used for Fig. 9 and Fig. 13.
For Fig. 9, which includes both type (iii)a and (iii)b simulations, we used the procedures described above, ex-cept we usedτ = 0.4, σ = 0 along the entire interface. by using σ f thres −f slip (a uniform distribution with support of extent f supp has standard deviation σ = f supp /(2 √ 3)) provides a data collapse under changes to k jun , f slip and f thres . Figure 18a shows unscaled data for uniform φ. To show that the scaling in equation (C2) is not limited to this particular shape of φ, we include data for bell-shaped φ in Fig. 18b. When we apply the scaling, in Fig. 18c, each distribution shape gets a data collapse. The master curves for the two distribution shapes are slightly different.
A point which we have used in Section V C is that wider φ correspond to lower effective local static friction. This can be seen directly in Fig. 18. Finally, we note that the scaling of µ eff s is the same as the scaling of prestress τ , a scaling that was also used in [9,10,76,77]. In this appendix we present the data underlying Fig. 11c. Figure 19 shows the full evolution of the events corresponding to a range of junction force distribution widths σ and withτ 0 = 0.05. For the lowest σ (strong interfaces) the front arrests, which is also reflected in Fig. 7. For the fronts that do reach the leading edge, the fronts are gradually changing from slow to fast propagation. These events plot in the "Slow" region of the front type map in Fig. 7 because they would arrest if the slow front mechanism was turned off. For each simulation we define the extent of the slow front region by visual inspection of the panels in Fig. 19. We then measure the front speed in the same way as usual, described in Appendix B. From these data we determine the average, minimum and maximum slow front speed values for these events. to 0.528 in steps of 0.059; the event for σ = 0 is not shown, but it too arrests early on. For each panel, the extent of the slow front was determined by visual inspection. On the left, the slow front starts after the step (ca. x = 0.04 m). On the right, because the transition to a fast front is more gradual, the boundary was set where the trend of near constant front velocity is broken. The exact limits chosen can be seen from the horizontal extent of the data in Fig. 11c. Panels (a-c) show a larger time interval than the rest of the panels. No data for Fig. 11c were extracted from the bottom right panel: even though the front has a slow part and plots in the slow part of the front type map, the slow part is too intermixed with the transition to fast rupture for the analysis to be carried out.