Synchronization of megathrust earthquakes to periodic slow slip events in a single-degree-of-freedom spring-slider model

Recently recognized slow slip events (SSEs) recurring in the deeper extensions of seismogenic zones along plate boundaries are drawing attention to their potential for triggering megathrust earthquakes that rupture the entire seismogenic zone. We describe how earthquakes simulated in a single-degree-of-freedom model are synchronized to the rhythm of imposed periodic SSEs. The time lag tQ from the one most recent SSE to a seismic event varies with system parameters and may take a broad range of eligible values between 0 and TSSE (SSE recurrence period). Earthquakes were found to synchronize with SSEs in various patterns depending on the proportion of SSE-driven loading within an SSE cycle, the recurrence period of the SSEs, and the duration of the SSEs, although synchronization itself remained a prevalent feature. Asynchrony was found only for long SSE durations.

www.nature.com/scientificreports www.nature.com/scientificreports/ evolution of the stress and the strength leading to synchronization. However, through experiments with various patterns of SSEs, we could find that the overall behavior of this synchronization phenomenon makes good sense in a lot of ways.
An earlier study pointed out that synchronization occurs in a spring-slider model 15 . It was numerically shown how more than one block engages in stick slip with one and the same period in a system of weakly interacting blocks that are mutually coupled via springs and have much the same characteristic periods when they are free from the influence of the other blocks. Another study attributed the observed clustering in time of great earthquakes in a given region to a synchronization phenomenon 16 . The situation assumed in the model study corresponds to that case. Unlike these studies, however, our study deals with earthquakes that are affected by a phenomenon (SSEs) with a much shorter period than the earthquakes themselves. We demonstrate that synchronization still occurs in that case.

Results
Model. To simulate plate-boundary megathrust earthquakes, we assume a system of a block being pulled on a frictional floor via a spring at velocity V load (Fig. 1a). Recurrent SSEs in the deeper extension of an earthquake fault are modeled as changes in the load-point velocity ( Fig. 1b; pulling displacement u load ; V load = du load /dt) so the impact of the repeating SSEs is integrated into the loading history. We conduct a quantitative study of dependence of the system behavior on the proportion r of the loading driven by SSEs within an SSE cycle, their recurrence period T SSE and their duration d SSE . An SSE begins at time t = jT SSE − d SSE and ends at t = jT SSE (j: integer). The V load is set at rV pl T SSE /d SSE when an SSE is going on and at (1 − r)V pl T SSE /(T SSE − d SSE ) during the rest of the time so the long-term pulling velocity remains at V pl , where V pl = 5 cm/yr and r = 0-1. A step displacement Δu SSE = rV pl T SSE is imposed, however, when d SSE = 0. The time evolution of the slip is calculated according to the Method. synchronization of earthquakes to repeating sses. In this subsection, we illustrate the example of a numerical test with r = 0.5, T SSE = 5 yr and d SSE = 0 to outline the synchronization phenomenon observed. The red curves in Fig. 2a,b show the time evolution of the slip rate V and stress τ, respectively, of the block in the case of L = 0.0618 m. The step-like growths in V and τ during the interseismic period are attributable to SSEs that occur every 5 years. In this case, earthquakes were found always to recur at an interval of T = 85 yr. A steady-loading model (V load = V pl ) with the same L value gives, by contrast, T 0 = 83.87 yr. The episodic-loading model gave a longer earthquake recurrence interval for an identical long-term loading rate. Figure 2 also shows, in green and blue, the results for larger L values of 0.0626 m and 0.0636 m, respectively. In both cases, T remained unchanged at 85 yr. As we have stated above, T 0 grows continuously in proportion to L in the steady-loading model. When L = 0.0636 m, T 0 = 86.31 yr, which means the recurrence interval in the episodic-loading model remains fixed and is now, contrary to the case of L = 0.0618 m, smaller than T 0 . And interestingly, the persistent value of T = 85 yr is exactly 17 times T SSE , which means the earthquakes are synchronized to the periodicity of the SSEs. This suggests a so-called entrainment phenomenon is at work. The fact that earthquakes see the rhythm of SSEs means that the timings of individual earthquakes is affected not only by the most recent SSE but also by the earlier sequence of SSEs. Figure 3b top shows the response of T over a broader range of L. Plotted in the panel are all the 200 or so values of T that emerged during the time interval t = 10,000 yr-30,000 yr in each of the simulation runs with different L values, which were sampled at intervals of 2 × 10 −4 m. The light-blue lines in the figure show, for reference, the values of T 0 (L) in the steady-loading model. The results came in two patterns: regime (i), where a unique T value is determined for a given L, and regime (ii), where more than one T value is ascribable to a given L.
The interval L = 0.0616 m-0.0636 m, for example, falls in regime (i). This interval is centered on a value of L, which we shall call L 17 , at which the natural characteristic period T 0 is exactly 17 times the SSE recurrence period (T 0 (L 17 ) = 17T SSE ). In this interval, which we shall call I 17 , T remains fixed at 85 yr = 17T SSE , and the periodicity of earthquakes is synchronized to the periodicity of SSEs. The cases shown in Fig. 2 all fall in the same interval I 17 . A range of greater L values, at L = 0.0638 m-0.0650 m, falls in regime (ii). More than one value of T, lying between 17T SSE and 18T SSE , emerges at each L during a single simulation run (we shall call this interval II 17−18 ). When L = 0.0644 m, for example, there are two alternating varieties of T, namely T 1 = 86.43 yr and T 2 = 88.57 yr, which add up to 35T SSE . In other words, an earthquake occurrence pattern is repeating with a period of P = T 1 + T 2 = 35T SSE , so earthquakes can be deemed synchronized to SSEs. The number of earthquakes per repeat cycle shall be called as n EQ . In this case, n EQ = 2; period-doubling bifurcation 17 occurred. When L = 0.0650 m, n EQ = 3 varieties of the recurrence interval, namely T 1 = 88.76 yr, T 2 = 89.55 yr and T 3 = 86.69 yr, repeat periodically. An earthquake occurrence pattern is repeating with a period of P = T 1 + T 2 + T 3 = 53T SSE . A relationship P = n SSE T SSE (n SSE : the number of SSEs per repeat cycle) always holds under regime (ii), so earthquakes can be deemed synchronized to SSEs under regime (ii), just as they are under regime (i). With L growing further, regime (i) sets in again at L = 0.0652 m-0.0672 m, with T fixed at 18T SSE (we shall call this interval I 18 ). This illustrates how regimes (i) and (ii) alternate with each other with changing system parameter L. With growing L, in other words, T continues to cling around T 0 (L) as it also grows, and goes alternately through an interval I m with T = mT SSE and an interval In both regimes (i) and (ii), P = n EQ T = n SSE T SSE , where T is the mean earthquake recurrence interval, so m = n SSE /n EQ , n EQ = 1 under regime (i) and n EQ > 1 under regime (ii). Regime (ii) is so-called high order synchronization (HOS) or n:m synchronization 18 . The triangles in Fig. 3b top show values of T , which grows monotonically, and in steps, with growing L. Apart from what looks like stair treads under regime (i)-a feature that is also seen in the variation pattern for T-the panel also shows the presence of fine, step-like features in T within regime (ii). This type of behavior is characteristic of synchronization phenomena and is known by the name of the devil's staircase 19 .
Let us introduce ΔT max , or the upper limit of |T k − T 0 | (amount of the deviation of an earthquake recurrence interval from the natural period T 0 ) within the entire experimented range L = 0.06-0.07 m as an indicator for measuring how far the values of T, realized under synchronization effects, can deviate from T 0 . The |T k − T 0 | was found to peak at both ends of each interval I m , and ΔT max = 1.51 yr = 0.30T SSE in Fig. 3b. The ΔT max can be interpreted as representing the magnitude of the "entrainment potential" of SSEs for entraining earthquake recurrence intervals into periodicity of a different nature.
The variation pattern for T may be interpreted as follows if we assume that SSEs cannot alter the periodicity of earthquakes above and beyond their entrainment potential. The T can only vary within the range T 0 www.nature.com/scientificreports www.nature.com/scientificreports/ − ΔT max < T < T 0 + ΔT max , and if any mT SSE values are available within that range, T is entrained into the closest of those mT SSE values, and regime (i) sets in. When no mT SSE value is available within that range, by contrast, the occurrence intervals keep varying within that eligible range and form a sequence of more than one seismic event, whose intervals sum up to an integer multiple of T SSE . That sequence serves as a repeat unit for achieving synchronization to SSEs. This is regime (ii). Regime (i) amounts, so to speak, to direct synchronization to the periodicity of SSEs, whereas regime (ii), with a weaker degree of synchronization than in (i), could be described as roundabout (indirect) synchronization. We here define a phenomenological parameter S, the proportion of the L intervals that realize regime (i), which is a measure of the prevalence of direct synchronization. In this case, S = 0.60.
Under the mechanism of entrainment proposed above, interval Interval I 18 , for example, covers the range L = 0.0652 m-0.0672 m. The widths of similar L ranges remained largely invariant regardless of the index m, and S = 0.60. That, combined with ΔT max = 1.51 yr = 0.30T SSE , confirms that Eq. (1) is met.
Let us next study the time lag t Q from the onset of the most recent SSE to a seismic event. In the case of Each of these timing relations between the SSEs and earthquakes remained invariant for all seismic events that recurred within a single simulation run. The earthquakes ended up phase-locked to the SSEs at the same T and t Q values even when the SSEs were shifted in occurrence time as part of the initial conditions. Figure 3b bottom shows t Q values over a broader range of L. As seen earlier in Fig. 2, t Q grows with increasing L within an interval I m along a slightly convex-downward curve shown in bold in the figure, and takes a broad variety of values ranging from 0 all the way up to T SSE . Within an interval II m−(m+1) , more than one T value corresponds to a given L, with each T k associated with its own, different t Q value. In both regimes (i) and (ii), the time lag t Q was found to take a broad range of eligible values between 0 and T SSE . The earthquakes are synchronized to the SSEs, without, however, showing any strong tendency to occur shortly following an SSE. Dependence of the system behavior on characteristics of the SSEs (their impact size r, recurrence period T SSE and duration d SSE ) will be studied in the following subsections.

Effect of r.
This subsection addresses dependence of the system behavior on the proportion r of the stress attributable to SSEs to all stress loaded on the fault during an SSE cycle. Apart from the above-described case of www.nature.com/scientificreports www.nature.com/scientificreports/ r = 0.5, we also studied the cases of r = 0.1 and 0.9, and the results are shown in Fig. 3a,c, respectively. The study found that earthquakes are synchronized with SSEs when r = 0.1 and 0.9 as well, with T and T generally growing with increasing L in following the "devil's staircase" pattern, except the case of r = 0.9, when regime (ii) does not exist. Note regime (ii) under r = 0.1 (Fig. 3a top) does exhibit small steps, though not seen with the coarse sampling of L in this figure (see Supplementary for finer sampling).
All recurrence intervals are close to T 0 when r is small. That is understandable, because a small r means only a small portion of the loads are driven by SSEs, so ΔT max , or the entrainment potential of the SSEs, is also small, and T, therefore, cannot deviate too far from the natural period T 0 .
It was also found that more T and t Q values correspond to a single L under regime (ii) when r = 0.1 than when r = 0.5, so n EQ and n SSE can sometimes take very large values. Among the largest, for example, are n EQ = 43 and n SSE = 719 for L = 0.0616 m. Phase-locking, however, still persisted in that case, and earthquakes were still found to be synchronized with SSEs. It appears likely that large n EQ values are required to form a repeat pattern under the constraint that, as stated above, T cannot deviate very far from T 0 when r is small.
When r = 0.1 and 0.9, t Q (Fig. 3 bottom) was found to grow with increasing L within an interval I m along a slightly convex-downward curve, just like in the case of r = 0.5. In both regimes (i) and (ii), the t Q values are distributed broadly across the range 0-T SSE , and the earthquakes exhibit no strong tendency to occur shortly following an SSE, even when SSE is large (r = 0.9).
Effect of T SSE . Let us next study how earthquake recurrence patterns vary with changes in the SSE recurrence period T SSE . We fixed r at 0.5 and d SSE at 0 as we changed T SSE from the 5 yr as before (Fig. 3b) to 1 yr and 10 yr. The results are shown in Fig. 4. Just like before, each T k was found to be associated with a fixed, unique value of t Q , and phase-locking was confirmed for all T SSE values. The T and T grow with increasing L as they go alternately through regimes (i) and (ii), and, just like in the earlier cases, the t Q values were found to be distributed broadly across the eligible range of 0-T SSE . It was found that S = 0.59, 0.60 and 0.59, and ΔT max = 0.30 yr, 1.51 yr and 3.14 yr for T SSE = 1 yr, 5 yr and 10 yr, respectively. In other words, ΔT max = 0.30T SSE , 0.30T SSE and 0.31T SSE , respectively, which confirms Eq. (1). www.nature.com/scientificreports www.nature.com/scientificreports/ The ratio ΔT max /T SSE remains almost unchanged, and that is also true in the cases of other r values, which we are not showing in our figures here. This, combined with ΔT max ∝ r found in the previous subsection, shows that ΔT max ∝ rT SSE . This means the entrainment potential is governed by the size of a single SSE step (rT SSE V pl ).
We further see that S does not depend on T SSE because T SSE dependences of numerator and denominator of Eq. (1) cancel out. A comparison between the cases of T SSE = 1 yr (Fig. 4a) and 5 yr (Fig. 3b) shows, for example, that the size of a single SSE step is five times smaller in the former case than in the latter, so the width of an L interval under regime I m is also five times smaller. By contrast, SSEs occur five times more often, which translates to five times more opportunities for entrainment into the SSE periodicity, so there are five times more I m intervals per unit L interval. Combined, S remains unchanged overall, and understandably, does not depend on T SSE .

Effect of d SSE .
We are finally studying changes in the SSE duration d SSE . With r fixed at 0.5 and T SSE at 5 yr, we studied the cases of d SSE = 1 yr and 2 yr along with the earlier case of d SSE = 0 (Fig. 3b). The results are shown in Fig. 5. The results for d SSE = 1 month or so, which are not shown in the figure, were hardly distinguishable from the results for d SSE = 0.
As readily seen from the existence of regime (i), synchronization occurred at least for certain L values even when d SSE > 0 (Fig. 5). However, in regime (ii), synchronization occurs only for some L values (black circles). For example, synchronization with n EQ = 7 is seen for L = 0.0684 m (black arrow) under d SSE = 2 yr (Fig. 5b bottom). For the other L values, synchronization does not seem to occur (red circles). For example, for L = 0.0612 m under d SSE = 2 yr (red arrow), no repeat pattern appeared throughout the entire experimented time (t = 10,000 yr-30,000 yr) worth 4,000 SSE cycles. We refer to this as asynchrony in the present study, although strictly speaking, we have not been able to make out for sure whether earthquakes and SSEs are still in sync in that case, only with immensely large n EQ and n SSE values beyond the tested number of SSE cycles, or they are really out of sync, with no periodic pattern whatsoever. For each L value with asynchrony, t Q took a broad range of eligible values.
It was found that S = 0.60, 0.38 and 0.14, and ΔT max = 1.51 yr (0.30T SSE ), 0.97 yr (0.19T SSE ) and 0.32 yr (0.06T SSE ), when d SSE = 0, 1 yr and 2 yr, respectively (Fig. 5 top). Equation (1) again holds true. The entrainment potential ΔT max decreases with increasing SSE duration, which, combined with findings from the previous subsection, shows that Figure 5 top www.nature.com/scientificreports www.nature.com/scientificreports/ shows the variation of T at each L is bounded by 2ΔT max constant throughout the experimented L range including both regime (i) and (ii), even when asynchrony occurs at some values of L in regime (ii). This suggests the concept of entrainment potential is valid even for the cases with asynchrony.
For regime (ii), we can define another measure of the prevalence of synchronization, S′; the proportion of the L values with synchrony among all the L values in the regime (ii). It was found S′ = 1.00, 0.45 and 0.09, when d SSE = 0, 1 yr and 2 yr, respectively. Asynchrony becomes more common with increasing d SSE .
Comparing Figure 5b (r = 0.5, T SSE = 5 yr and d SSE = 2 yr) and Figure 3a (r = 0.1, T SSE = 5 yr and d SSE = 0), we notice the entrainment potential ΔT max and S are about the same. Values of L with asynchrony, however, emerged only in the former case (S′ = 0.09), so it can be said the system is less prone to synchronization in the former case than in the latter. This suggests the presence of a certain effect of d SSE , on the prevalence of indirect synchronization in regime (ii).
In summary, the prevalence of direct synchronization (regime (i)), measured with S, is controlled by the entrainment potential ΔT max = f(d SSE )rT SSE . Additionally, the prevalence of indirect synchronization within regime (ii), measured with S′, decreases with increasing d SSE .

Discussion and Conclusion
We have used a spring-slider model to simulate the occurrence of earthquakes that are synchronized to the periodicity of loading by recurring SSEs. All seismic events in the experimented cases were found to be synchronized with SSEs when d SSE = 0, and they were also found always to be in sync even when the single SSE size was small, such as when r = 0.1 and T SSE = 5 yr. We could confirm synchronization down to r as small as 0.03 when T SSE = 5 yr (figure not shown). When d SSE > 0, by contrast, earthquakes were found to be out of sync with SSEs at certain L values. It is understandable that, the longer the SSE duration, the smaller the difference in the loading rate between when an SSE is going on and during the rest of the time, and hence the system becomes less prone to synchronization because an SSE becomes less of a special event or a clear time marker. When r = 0.5, d SSE = 2 yr and T SSE = 5 yr, for example, there is an ongoing SSE over 40% of all time, and the loading rate then is only 1.5 times the corresponding rate during the rest of the time, and asynchrony was found for the majority of the L values studied. Even in that case, however, synchronization was still found for some L values.
Two modes of synchronization were recognized when earthquakes were in sync with SSEs: regime (i) of direct entrainment into the periodicity of the SSEs, and regime (ii) of indirect entrainment, which involves the formation of a repeat pattern consisting of more than one seismic event. Periodic SSEs have the potential to shift an earthquake recurrence interval from the natural period T 0 by up to the magnitude of their entrainment potential ΔT max . Whether the regime will be (i) or (ii) is determined by whether any mT SSE values are available within that allowable range. Then ΔT max controls the prevalence of direct synchronization, whereas the prevalence of indirect synchronization in regime (ii) is affected by d SSE .
The settings we have studied in the present article-(1) a smaller r, (2) a smaller T SSE and (3) a larger d SSE -can all be regarded as a move toward the limit of the steady-loading model. The earthquake recurrence intervals T approached T 0 under settings (1) and (3), and the system behavior approached that of the steady-loading model, wherein earthquake recurrence intervals vary continuously and linearly with L. No such approach, however, occurred under setting (2), because S does not depend directly on T SSE . Case (2) highlights a nonlinear nature of the synchronization system being studied here.
Our simulations over a wide range of settings have robustly exhibited synchronization of earthquakes to periodic SSEs. The appearance of synchronization means that the whole history of SSEs matters for the timing of individual earthquakes. Time delay of an earthquake from the most recent SSE is distributed broadly across the entire range of eligible values 0-T SSE . Thus, from the viewpoint of forecasting, the risk enhancement implied by the observation of an SSE is much less than that expected from the clock advance effect due to a single SSE imposed at an arbitrary timing 12 .
Future research could provide evidence of synchronization between earthquakes and SSEs in field data. The present article has shown that earthquakes may occur not necessarily at fixed recurrence intervals (regime (i)) but also in a complicated pattern of synchronization that involves more than one recurrence interval, so that may not give the impression at first glance that the seismic events are synchronized with SSEs (regime (ii)). As long as few SSEs and earthquakes have been observed, it remains difficult to discuss the relations between the SSEs and seismic events unless a clear spatiotemporal relationship is recognizable as in the Guerrero case mentioned in the introduction 10 . In the future, however, the complicated synchronization phenomena may be found in the modern observation data accumulated over a long time.

Method
Under a pulling displacement u load (Fig. 1b), the equilibrium of forces acting on the block (Fig. 1a) can be described quasi-dynamically with the following Eq. (2) of motion 20,21 by using the spring constant k, the block displacement u, the block slip rate V and the friction τ: where the S-wave velocity V s = 3.27 km/s and the rigidity G = 30 GPa. We assume the friction τ follows a rate-and-state friction (RSF) law (Eq. (3)) and use the aging law (Eq. (4)) to serve as an evolution law for the strength parameter Φ 22 : n www.nature.com/scientificreports www.nature.com/scientificreports/ where a, b and L are friction parameters, σ n is the effective normal stress on the block's sliding face, and V * is an arbitrary reference velocity (we set V * = V pl in the present study). The τ * refers to the frictional stress under the steady state with velocity V * . By way of initial condition, we set V at 0.9V pl and Φ at the corresponding steady-state value when t = 0. The first and second terms on the right-hand side of Eq. (4) represent the processes of time-dependent healing and slip weakening of the strength, respectively. The block motion in this system is generally takes the form of intermittent stick slip (earthquakes) when a < b and k < k critical (≡ σ n (b − a)/L). We set σ n at 100 MPa; a = 0.010 and b = 0.012; and k at 0.5k critical , because we are interested in earthquake occurrences. We use the Runge-Kutta method with variable time steps 23 to solve the system of Eqs (2)(3)(4) to obtain time evolutions of the slip rate and the stress.

Data Availability
No datasets were used in the current study.