Aftershocks in slowly compressed bulk metallic glasses: Experiments and theory.

We observe two distinct interevent time patterns in the slip avalanches of compressed bulk metallic glasses (BMGs). Small slip avalanches cluster together in time, but large slip avalanches recur roughly periodically. We compare the timing patterns of BMG slip avalanches with timing patterns of earthquakes and with the predictions of a mean-field model. The time clustering of small avalanches is similar to the known time clustering of earthquake foreshocks and aftershocks.


I. INTRODUCTION
Bulk metallic glasses (BMGs) are amorphous alloys with superior mechanical properties compared to many conventional structural materials; however, they typically break in abrupt brittle failure [1][2][3]. Slowly compressed specimens of BMGs deform elastically to high stresses, but when deformation enters the plastic regime, BMGs deform in a jerky manner with sudden stress drops ("slips" or "avalanches"), culminating in catastrophic fracture [4][5][6][7][8][9][10][11][12]. We extract interevent time patterns between these stress drops that demonstrate the occurrence of quasiperiodic large slips, time-clustered small slips, and foreshocks and aftershocks of both large slips and small slips in BMGs. The observed BMG slip statistics are similar to earthquake statistics, and small slip avalanches even show aftershocks and foreshocks akin to earthquakes. The interevent time patterns from experiments are compared to the predictions of a simple mean-field model. The results are important for materials testing and hazard prevention.
The equipment and procedure used to gather the data in this experiment have been detailed previously [7,13]. Two mmscale BMG specimens (Zr 45 Hf 12 Nb 5 Cu 15.4 Ni 12.6 Al 10 ) were loaded into a mechanical test system that increased the compressive strain on each at a constant displacement rate to achieve a nominal strain rate of 10 −4 s −1 [7]. The stress was recorded with a time resolution of 10 μs. The stress vs time traces were Wiener filtered to reduce high-frequency noise, and the start times (end times) of slip avalanches were identified with local maxima (minima) in the Wiener-filtered stress vs time traces [7]. This avalanche-detection algorithm uses a velocity threshold of 0 MPa/s. Previous work has * Retired. † dahmen@illinois.edu shown that there are two regimes of slip avalanches: small avalanches obey power-law scaling, while large avalanches exhibit many differences from small avalanches, for example, in the typical avalanche shape [7]. We estimate 10 MPa to be the approximate classification boundary between small slip avalanches and large slip avalanches. Specimen (1) yielded 268 small avalanches (0.2 MPa avalanche size < 10 MPa) and 63 large avalanches (10 MPa avalanche size). Specimen (2) yielded 338 small avalanches and 53 large avalanches in the same small and large size ranges. We analyzed the slip avalanches from both BMG specimens, and we then aggregated the results into a single dataset in this paper.
The slip statistics of the filtered time-dependent stress were previously shown to agree with 12 different statistical predictions of a simple mean-field model [7]. In addition, as predicted by the model, two separate avalanche regimes were apparent both in the complementary cumulative distribution function (CCDF) of avalanche sizes and in the dynamics of individual avalanches. Small, localized avalanches propagate intermittently in a jerky manner. In contrast, large, specimen-spanning avalanches have long incubation periods, with smooth and sharply peaked stress rate vs time profiles for individual events [7,14]. In this work, we show that while the large avalanches recur quasiperiodically, the small avalanches cluster together at short timescales and occur randomly at long timescales.

II. INTEREVENT TIME CORRELATIONS
The temporal correlations between different avalanches contain important information about the plastic deformation process of BMGs [15]. We start by examining the probability distribution of waiting times for avalanches in the smalland large-avalanche regimes. Interevent times are defined as the difference between nucleation times, or "start times," of temporally neighboring avalanches. Figure 1(a) shows the probability density of interevent times computed from the set of small avalanches, with 0.2 MPa small avalanche size < 10 MPa. The resulting probability density of interevent times does not follow a simple functional form, but the small-avalanche interevent times are heavily weighted toward times less than FIG. 1. Interevent times for slip avalanches experimentally measured in mm-scale specimens of bulk metallic glass slowly compressed at a strain rate of 10 −4 s −1 . (a) Probability density of interevent times between small slip avalanches. The probability density between 2 and 30 ms is consistent with a power law of exponent −1.4 ± 0.2. (b) Histogram of interevent times between large slip avalanches. 1 s. Figure 1(b) shows the histogram of interevent times computed from the set of large avalanches with sizes larger than 10 MPa. In contrast to the small avalanches, the large-avalanche interevent time distribution has its maximum around 7 s and another local maximum around 3 s, indicating a quasiperiodic pattern. Appendix B shows that the two BMG specimens have different characteristic interevent times, with the result that the combined dataset displays the bimodal structure observed in Fig. 1(b). The characteristic interevent times of the two specimens may depend, e.g., on different distributions or locations of shear bands between the two specimens.
To test for precursors to large slips, for temporal clustering, and for temporal quasiperiodicity, we perform the Bi test on the avalanche start times [16][17][18]. The Bi test was originally developed to detect correlations in absorption lines in quasars, and it has since been applied to find correlations in avalanche processes. A dataset is said to "pass" the Bi test if the dataset is consistent with a Poisson process. If the dataset fails the Bi test, then the data are inferred to be inconsistent with a Poisson process. The Bi test has two powerful statistical advantages in detecting correlations. First, the Bi test is robust to both homogeneous and inhomogeneous Poisson processes-either type of process will pass the Bi test. Second, if a dataset fails the Bi test, then the way in which the dataset fails the Bi test yields important insights into the temporal correlations of the underlying process [18]. Datasets that exhibit clustering fail the Bi test in a qualitatively different way than datasets that exhibit quasiperiodicity.
The Bi test applied here takes the chronological array of avalanche start times as input. Each avalanche, indexed by integer k, from the third to the third-to-last avalanche, has a value H k associated with it as defined below. For each avalanche, δt k and δτ k are defined as where t k is the start time of the kth avalanche. For the kth avalanche, H k is then defined by For a Poisson process, either homogeneous or inhomogeneous, the cumulative distribution function (CDF) of the set of H k values fluctuates around a straight line from (0,0) to (1,1) [16]. Datasets that exhibit quasiperiodicity have a sudden increase in the CDF of H k values at H k = 2/3. In the limit of perfect periodicity (for which δt k = δτ k ), all H k = 2/3. In contrast, datasets that exhibit clustering have anomalously short times between avalanches within a cluster and anomalously long times between avalanches in separate clusters. The corresponding CDF in this case exhibits a rapid increase near small and large H k values. The Bi test therefore has a qualitative interpretation yielding a quasiperiodic process, a Poisson process, or a clustering process, as indicated by the shape of the CDF curve [18]. This interpretation becomes quantitative with the inclusion of a two-sided Kolmogorov-Smirnov (KS) test [17]. This two-sided KS test takes two input values: the maximum absolute vertical deviation of the CDF of H k values from the hypothetical straight line, and the number of H k values. The output of the KS test is the probability that a Poisson process with the given number of data points would yield such a CDF, with such a maximum deviation from the hypothetical straight-line CDF. The null hypothesis of this test is that the process is Poissonian. The Bi test cannot prove that a process is Poissonian, only that a process is statistically consistent with a Poisson process. Conversely, the Bi test can show that an underlying process is not Poissonian to a high degree of probability.
The results of applying the Bi test to small and large avalanches are shown in Bi test for small (blue ×) and large (red circle) avalanche interevent times experimentally measured in slowly compressed bulk metallic glass. The solid line represents the ideal Bi test for a local Poisson process. The curve for the small-avalanche Bi test rapidly increases near small and large values of H k ; this shape is typical of correlated times that occur in clusters. The curve for the large-avalanche Bi test rapidly increases near H k = 2/3; this shape is typical of correlated times that occur at fairly regular intervals. Small avalanches cluster together; large avalanches recur nearly periodically. P = 1.1 × 10 −12 of being Poissonian (effectively zero in both cases). Moreover, the two avalanche regimes fail the Bi test in strikingly different ways. The Bi test for large avalanches exhibits a rapid increase at approximately H k = 2/3 on the horizontal axis; this behavior is characteristic of a quasiperiodic process. The Bi test for small avalanches instead exhibits a rapid increase near small and large H k values; this behavior is characteristic of a clustering process. The Bi test therefore indicates that large slip avalanches occur approximately periodically, while small slip avalanches tend to cluster together in time. In Appendix D, we also apply the Bi test to the combined set of all avalanches, both large and small. Since most of our avalanches are small avalanches, the combined Bi test is similar to the Bi test for small avalanches only.

III. AFTERSHOCK AND FORESHOCK RATES
To determine the timescale of the small-avalanche clustering, we borrow the concepts of "aftershocks" and "foreshocks" from earthquakes and apply these concepts to BMG slip avalanches. We define the aftershock sequence for a given avalanche (known as a mainshock) to consist of all slip avalanches that occur after the mainshock, until an avalanche with size greater than the mainshock is reached. For every avalanche and its associated aftershock sequence, we define the time origin (t = 0) to be the start time of the mainshock. The start time of each aftershock is then the time since the beginning of the associated mainshock. These definitions of aftershocks and aftershock start times are consistent with prior work [17,[19][20][21]. We define the aftershock rate to be the number of aftershock nucleations per unit time.
To calculate aftershock rates, we divide the time axis into logarithmically spaced bins. For each time bin, we count the number of aftershocks that start in a given bin, iterating over all aftershock sequences. We also count the amount of time in each bin during which an avalanche is not occurring and refer to this as the "dead time" in the bin. The number of avalanches that start in each bin divided by the dead time in each bin yields the aftershock rate in each bin. This definition arises from consideration of the sample space of detectable start times. Our algorithm cannot detect an aftershock that nucleates while another aftershock is still ongoing; i.e., there is no opportunity to detect or record the start time of the second avalanche, since the two (or more) avalanches merge into a single large avalanche [22]. In Appendix E, we explore the effects of ignoring this issue by using the "real time" rather than only the dead time. The "real time" is the full amount of time that each aftershock sequence spends in each time bin, regardless of whether an avalanche is occurring or not. This definition of aftershock sequences does not strictly delineate avalanches as either mainshocks or aftershocks; strict classifications cannot be made since we do not measure the spatial locations of avalanches within the BMG specimen [23]. Such a classification for earthquakes uses a combination of temporal and spatial distance between earthquake events; however, our measurements do not include spatial information. We instead investigate the average avalanche rate following an initial avalanche for the series of subsequent avalanches that are smaller than the initial avalanche.
We find that both small and large avalanches have aftershocks. Figure 3 shows the average rate of aftershocks after FIG. 3. Average aftershock rates for small (blue ×) and large (red circle) mainshocks experimentally measured in slowly compressed bulk metallic glass. The power-law region of the aftershock rate of small mainshocks is consistent with a power-law exponent of −1.3 ± 0.15. The time origin is the beginning of the mainshock. After a small mainshock, the aftershock rate follows a power law as the rate decays with time; the deviation from power-law behavior at short times is likely due to the minimum detectable avalanche size. After a large mainshock, the aftershock rate follows a steep curve as the rate decays with time. The two points that have no error bars were calculated from bins in which there was only one detected aftershock nucleation and should be thought of as having relatively large uncertainty. small avalanches and after large avalanches. The aftershock rate following a small avalanche exhibits a much smoother power-law decay than the aftershock rate following a large avalanche. The aftershock rate after small avalanches decays as a power law with exponent −1.3 ± 0.15, ending before 100 ms since the mainshock. A power-law dependence of the aftershock rate vs time, with an exponent close to 1, has also been observed in earthquakes. The Omori law describes a similar aftershock rate vs time in earthquake statistics [24].
Foreshocks and foreshock sequences are defined in a similar way as aftershocks and aftershock sequences, but with avalanches that occur prior to the given avalanche (mainshock) instead of after the given avalanche. The time origin of each foreshock sequence is the start time of the associated mainshock; the foreshock sequence then progresses toward earlier times and toward earlier avalanches, and the foreshock sequence ends when it reaches an avalanche of size larger than the associated mainshock. The foreshock rate is the number of foreshock nucleations per unit time. Figure 4 shows the average foreshock rate prior to small avalanches and prior to large avalanches. In contrast to the aftershock rates, the foreshock rates for small and large avalanches are similar. The foreshock rates for both small and large avalanches roughly follow a power law with exponent −1.15 ± 0.15, valid only within <100 ms of the mainshock. The power-law growth of the foreshock rates in BMG slip avalanches is also similar to power-law growth of the foreshock rates in some earthquake catalogs [24].

IV. DISCUSSION AND CONCLUSION
The interevent time patterns of slip avalanches in BMGs are consistent with the interevent time patterns of a simple mean-field model of plastic deformation [25,26]. Previous work shows that slip avalanches in BMGs are similar to slip avalanches in the mean-field model, across a dozen statistical quantities [7]. The mean-field model describes the dynamics of many shear-stressed cells, which are elastically coupled both to each other (through a mean-field elastic interaction) and to an externally imposed, slowly increasing shear stress. Each cell remains stationary until the stress on that cell exceeds a failure stress, causing that cell to move to a new position until the stress on that cell falls below an arrest stress. When a cell slips in this way, the stress on all the other cells increases toward the failure stress, via the mean-field elastic coupling between each pair of cells. A moving cell can therefore cause other cells to fail, which in turn cause more cells to fail, creating a chain-reaction avalanche process. When weakening is included in this model, whereby cells that have already moved are "weakened" and are easier to move again, the resulting dynamics include quasiperiodic large avalanches. These quasiperiodic large slips predicted by the mean-field model are evident in the data as discussed earlier. Additionally, aftershocks and avalanche clustering are present in an extension of the simple mean-field model, via inclusion of a slow relaxation time intrinsic to the material [27]. We note that a version of the mean-field model employed here was originally used to describe earthquakes [25,27,28].
We have shown that both small and large slip avalanches have aftershock and foreshocks. Small slip avalanches have power-law aftershock and foreshock rates that vary with time; power-law aftershock and foreshock rates are also present in earthquake sequences, as described by the Omori law for earthquakes. These experimental earthquakelike results complement prior simulations of glassy materials [29]. These results using avalanche start times also build on the previous analysis of a stressed Al3%-Mg alloy, in which the existence of temporal correlations of acoustic emission (AE) events can be inferred through multifractal analysis of the coarse-grained AE signal [30]. In addition, prior work has shown that the smallavalanche energy distribution follows a power law, which is similar to the Gutenberg-Richter scaling law for earthquake magnitudes [14]. Furthermore, BMG slip avalanches exhibit these interevent time patterns across more than two orders of magnitude in time: small avalanches display clustering on timescales down to a few milliseconds (Fig. 3), while large avalanches display periodicity on a timescale of approximately 7 s. The presence of time-series patterns across different temporal orders of magnitude is akin to the wide range of time scales for which temporal correlations are present in earthquake catalogs [31,32].
One limitation of our measurement is that we cannot detect the smallest avalanches; avalanches of size less than 0.2 MPa are lost in the noise of our measurement. Our conclusions for small avalanches are therefore valid only for the "restricted" small-avalanche regime, i.e., for small avalanches with 0.2 MPa avalanche size < 10 MPa. One might ask if the inclusion of the smallest avalanches below the noise floor may change the observed small-avalanche temporal clustering demonstrated here. We note, however, that even if inclusion of the smallest avalanches would yield an uncorrelated Poisson process, the observed "restricted" small-avalanche temporal clustering shown here is important both theoretically and practically. In real systems, for example, earthquakes, avalanches  that are larger than some noise floor that exists in most real systems are of more practical consequence. Even if the smallest slip avalanches yield different behavior, we have found that the more important, relatively larger scaling avalanches are clustered together in time. We have also tested the effect of dropping the smallest avalanches on the test in our mean-field simulations and found no significant effect [33].
A recent paper discussed the possibility that observed temporal correlations in an avalanche process could be spurious artifacts of the avalanche-detection algorithm, arising from the use of a finite stress-rate threshold [34]. Since we use Wiener filtering to greatly reduce the high-frequency noise, we are able to use a stress-rate threshold of 0 MPa/s to detect possible avalanches. Only the avalanches for which the total stress change is at least 0.2 MPa are kept [7]. Our avalanche-detection algorithm thereby avoids this finite-rate-threshold problem. Nevertheless, we found the analysis in [34] to be useful and interesting. In Appendix F, we use a similar idea inspired by that paper to analyze our data, and we reach the same conclusions as presented here.
Additionally, the aforementioned paper shows that the threshold-induced spurious correlation effect results in a particular power-law exponent of the interevent time distribution. Following the notation of [34], if we denote the power-law exponent of the interevent time distribution by τ Tw , and if we denote the power-law exponent of the avalanche duration distribution by τ T , then according to [34] the presence of spurious correlation effects would imply that τ Tw = τ T . In our previous work, Fig. S8 in the Supplemental Material shows that the slip avalanches in our BMG experiment exhibit a power-law exponent τ T = 2 in their duration distribution (that figure shows that the CCDF of avalanche durations follows a power law of −1, and so the probability distribution of avalanche durations follows a power law of −τ T = − 2) [7]. Therefore, threshold-induced spurious correlations would yield a power-law regime of avalanche interevent times with scaling exponent τ Tw = 2. This exponent is not present in our interevent time distribution shown in Fig. 1(a), and so our results are incompatible with a threshold-induced spurious correlation.
This study shows that the experimental data on slowly compressed BMGs not only agree with the slip statistics predicted by a simple mean-field model, as shown previously, but also with important time-series properties related to aftershocks and foreshocks of earthquakes. The mean-field model has been applied to a broad array of deformation phenomena, across length scales ranging from the geologic to the microscopic [25][26][27]35]. In the future, it would be interesting to study the aftershock and foreshock properties in these other systems as well in order to test the generality of our findings. The results presented here already suggest an intriguing link between the interevent times of BMG slip avalanches and those of earthquakes that warrants further consideration of BMGs as model laboratory systems for studying earthquakes. In this paper we use 10 MPa as the classification boundary separating small avalanches from large avalanches. This classification boundary is not strict; it is possible that some small "scaling" avalanches exceed this size, and it is possible that a few large system-spanning avalanches fall below this size. Nevertheless, this simple classification boundary yields a robust difference in the temporal correlations. In the Supplemental Material of our previous work [7], the size CCDF is shown for each BMG specimen separately, and a size boundary that separates large (system-spanning) avalanches from small (scaling) avalanches can be identified as roughly 10 MPa. Figure 1(b) shows the interevent times between large avalanches (avalanches with stress drop size >10 MPa) for the full set of interevent times aggregated from both specimens. There is an apparent bimodal distribution in this figure, with one peak near 3 s and another peak near 7 s. We investigate the possibility that this double peak reflects different characteristic interevent times for the two different specimens. Moreover, visual inspection of Fig. 5(a) suggests that the characteristic interevent time increases during the second half of the time trace. We investigate the possibility that the avalanche interevent time distribution evolves with time. Figures 6(a) and 6(b) show the probability density of interevent times between small avalanches, separated both by specimen and by temporal half of the stress vs time trace. These figures are to be compared with the aggregated data shown in Fig. 1(a). The probability density of small-avalanche interevent times is stable and shows minimal changes as a function of the specimen or temporal half of the stress vs time trace.

APPENDIX B: SPECIMEN DEPENDENCE AND TIME DEPENDENCE OF THE DISTRIBUTION OF AVALANCHE INTEREVENT TIMES
FIG. 6. Probability density of interevent times between small slip avalanches for (a) BMG specimen (1) and (b) BMG specimen (2). The curves labeled "first half" ("second half") are only calculated from the interevent times from the temporal first half (second half) of the stress vs time trace. The small-avalanche interevent time distribution is relatively constant across time and from one specimen to the other. The few points in (a) and in subsequent figures that have no error bars were calculated from bins in which there was only one interevent time and should be thought of as having relatively large uncertainty. Figures 7(a) and 7(b) show the histogram of interevent times between large avalanches, separated both by specimen and by temporal half of the stress vs time trace. These figures are to be compared with Fig. 1(b). The histogram of large-avalanche interevent times in Fig. 7(a) shows that the characteristic interevent time changes from about 3 s in the first half of the specimen (1) experiment to about 6 s in the second half, while the characteristic interevent time is stable at 7 s in specimen (2). The origin of the discrepancy between the large-avalanche interevent times is not determined here, but these large events are thought to be system-spanning avalanches that depend on the boundaries of the BMG specimens [7]. Both of these characteristic timescales, 3 and 7 s, are orders of magnitude longer than the small-avalanche interevent times shown in Fig. 1(a).  (1) and (b) BMG specimen (2). The curves labeled "first half" ("second half") are only calculated from the interevent times from the temporal first half (second half) of the stress vs time trace. The characteristic large-avalanche interevent time changes by about a factor of 2 from the first temporal half of specimen (1) to the second temporal half. In contrast, the characteristic interevent time is more stable for specimen (2).

APPENDIX C: SPECIMEN DEPENDENCE AND TIME DEPENDENCE OF THE Bi TEST RESULTS
As in the preceding section, we investigate whether the results of the Bi test depend on the specific BMG specimen and also whether the results differ between the first half and the second half of each specimen's stress vs time trace. Figures 8(a)  and 8(b) show the Bi tests for large and small avalanches, separated by specimen, and separated by temporal half of the stress vs time trace. These figures are to be compared with Fig. 2 in the main text. Despite the observed change in characteristic interevent time in specimen (1) [shown in Fig. 7(a)], the main results of the Bi test (i.e. small-avalanche clustering and large-avalanche quasiperiodicity) are the same for both specimens and do not appreciably change with time.

APPENDIX D: Bi TEST APPLIED TO THE SET OF BOTH LARGE AND SMALL AVALANCHES
In the main text, we applied the Bi test separately to the set of interevent times between small avalanches and between large avalanches. In Fig. 9 we show the Bi test applied to 063005-7 FIG. 9. Bi test applied to the set of all detected avalanches using all interevent times aggregated from both temporal halves of both specimens. These interevent times are not separated either by temporal half of the experiment or by avalanche size into small vs large avalanches as in the main text; instead these interevent times are between any avalanches with avalanche size above the noise floor (0.2 MPa avalanche size). the set of interevent times between the full set of avalanches, both large and small. Since most of the avalanches are small avalanches, as shown in our previous publication [7], we expect the resulting all-avalanche Bi test to be comparable to the small-avalanche Bi test in the main text. Indeed, the all-avalanche Bi test indicates a clustering pattern similar to the small-avalanche clustering in Fig. 2.

APPENDIX E: AFTERSHOCK AND FORESHOCK RATES CALCULATED WITHOUT ONLY USING THE "DEAD TIME"
In the main text, we explain how we calculate aftershock rates and foreshock rates using the "dead time," i.e., the time during which an avalanche is not occurring. Our data are not spatially resolved; we measure stress vs time effectively averaged across the entire specimen. Therefore, if a separate avalanche starts while a previous avalanche is still ongoing, the two avalanches merge to appear as one avalanche in our data. The result is that any time during which an avalanche is occurring is not in the "sample space" of detectable start times. In calculating the aftershock or foreshock rates, we divide the number of avalanche start times by the time elapsed during which we are able to detect avalanche start times. The time during which we are able to detect avalanche start times is the "dead time." For the sake of comparison, we investigate an alternative definition of aftershock and foreshock rates. With this alternative definition, we calculate the rates by using the "full time" of each aftershock sequence instead of only the (shorter) dead time; the "full time" is the total amount of time that each aftershock sequence spends in each time bin regardless of whether or not an avalanche is occurring. The results are shown in Figs. 10(a) and 10(b). With this alternative definition, the power-law behavior of aftershock rates and foreshock rates of small avalanches is similar to the power law of Fig. 3 in the main text, within the exponent error bars. With this analysis, the power-law exponent of the aftershock rate after a small avalanche is −1.2 ± 0.15, and the power-law exponent of the foreshock rate before a small avalanche is −1.15 ± 0.2. With this alternative definition, aftershock and foreshock rates still follow Omori-like power laws as a function of time, so that the analogy to earthquake behavior is still observed.

APPENDIX F: COMBINING AVALANCHES CLOSE TOGETHER IN TIME
As mentioned in the main text, a recent paper explored the pitfalls of extracting temporal correlations from avalanche data that were analyzed using a nonzero-velocity threshold [34].  [7]. The shape collapse quality becomes worse as the minimum time increases.
Our data have been Wiener filtered, so that the high-frequency noise has been strongly attenuated. Wiener filtering enables us to use a zero-velocity threshold to detect avalanches, thereby preventing the described problem.
Nevertheless, in the spirit of that interesting paper, we consider a related process by which spurious temporal correlations could appear. Additive noise on top of the "real" avalanche velocity signal could cause the avalanche velocity to appear 063005-9 FIG. 12. Avalanches analyzed with a "minimum time" of 10 −3 s between the ending time of any avalanche and the beginning time of the subsequent avalanche. (a-e) show the same quantities that are reported in the main text, but with a nonzero "minimum time" for these figures. The general conclusions of the main paper, i.e., that small avalanches cluster together while large avalanches recur quasiperiodically, hold true even with this modified avalanche analysis.
FIG. 13. The Bi test using a delayed Poisson process as the null hypothesis [36]. (a) The delayed Poisson process Bi test applied to the same data as in Fig. 2 (large vs small avalanches). (b) The delayed Poisson process Bi test applied to the same data as in Fig. 8(a). (c) The delayed Poisson process Bi test applied to the same data as in Fig. 8(b).
to drop below 0 MPa/s for a small duration while a single avalanche is still ongoing. Such a downward noise spike could therefore cause a single avalanche to incorrectly appear to be multiple closely spaced avalanches, thereby leading to spurious temporal clustering of avalanches.
We devise the following method to explore whether this effect causes the observed pattern of small-avalanche clustering. We reanalyze the filtered stress vs time traces in order to detect avalanches, but this time we introduce an adjustable quantity that we call the "minimum time" between avalanches. The "minimum time" is the minimum amount of time allowed between the ending time of any one avalanche and the beginning time of the subsequent avalanche. If two avalanches are thus too close together in time, we merge them into one avalanche by keeping only the earlier beginning time and the later ending time. In this way, we recombine avalanches that have possibly been erroneously broken up into multiple avalanches; however, if the "minimum time" is too large, then we may instead mistakenly combine avalanches that truly are separate avalanche events and therefore destroy the real temporal clustering.
We perform this analysis for a range of "minimum times," extending from 0 s (which is identical to the analysis in the main text) to 0.1 s. When the value of this "minimum time" quantity becomes too large, the avalanche shapes will be distorted due to the inclusion of too many lengths of near-zero velocity. For each value of the "minimum time," we sort the avalanches into bins based on size, and then take the average avalanche shape in each bin. We then rescale the average avalanche shape in each bin by the expected mean-field size-dependent scaling factor, and we compare these size-scaled shapes to the expected sizebinned mean-field avalanche shape in Figs. 11(a)-11(f). The largest "minimum time" that exhibits the expected mean-field shape collapse is shown in Fig. 11(d), with a minimum time of 10 −3 s. Beyond that value, the average avalanche shapes are diminished and distorted due to the inclusion of long stretches of near-zero velocity. We conclude that the largest plausible value of the "minimum time" is 10 −3 s, and that for larger values we are merging separate avalanche events that should not be combined into one event.
Using this maximum "minimum time" value, wherein we combine avalanches that are separated in time by less than 10 −3 s, we recalculate the figures from the main text to see if the prior observed patterns still appear with this modified analysis [Figs. 12(a)-12(e)]. The temporal correlations are still present but are slightly weaker on timescales close to the "minimum time," as is expected when the number of interevent times close to the "minimum time" is reduced by combining multiple avalanches into a single avalanche event. We again find that small avalanches cluster together in time, while large avalanches occur quasiperiodically. Small avalanches continue to exhibit Omori-like aftershock and foreshock power-law behavior. We conclude that the type of spurious noise-spike effect described above does not account for the range of earthquakelike temporal patterns that we observe in BMG slip avalanches.

APPENDIX G: Bi TEST WITH A DELAYED POISSON PROCESS NULL HYPOTHESIS
The Bi test as described above uses a Poisson point process as the null hypothesis. If a process is not pointlike, i.e., if the avalanche durations are not small compared to the interevent 063005-11 times, then the Bi test may be adapted to use a delayed Poisson process as the null hypothesis [36]. In a delayed Poisson process Bi test, each interevent time is calculated as an avalanche start time minus the previous avalanche end time. This is in contrast to the standard Poisson process Bi test, in which each interevent time is calculated as an avalanche start time minus the previous avalanche start time. These two definitions are equivalent for processes that have only zero-duration avalanche events. The delayed Poisson Bi test and the standard Bi test yield very similar results for the BMG avalanches analyzed in this paper. Figure 13(a) shows the delayed Poisson Bi test applied to the same data as in Fig. 2. The small avalanches have a probability P = 1.9 × 10 −52 of being Poissonian, and the large avalanches have a probability P = 1.1 × 10 −12 of being Poissonian. Figure 13(b) shows the delayed Poisson Bi test applied to the same data as in Fig. 8(a). The small avalanches from the first half (second half) of specimen (1) have a probability P = 9.5 × 10 −12 (P = 7.3 × 10 −20 ) of being Poissonian. The large avalanches from the first half (second half) of specimen (1) have a probability P = 5.5 × 10 −6 (P = 0.014) of being Poissonian. Figure 13(c) shows the delayed Poisson Bi test applied to the same data as in Fig. 8(b). The small avalanches from the first half (second half) of specimen (2) have a probability P = 6.0 × 10 −16 (P = 1.9 × 10 −9 ) of being Poissonian. The large avalanches from the first half (second half) of specimen (2) have a probability P = 2.4 × 10 −5 (P = 0.0018) of being Poissonian.