Synchronization of coupled stick-slip oscillators

A rationale is provided for the emergence of synchronization in a system of coupled oscillators in a stick-slip motion. The single oscillator has a limit cycle in a region of the state space for each parameter set beyond the supercritical Hopf bifurcation. The two-oscillator system that has similar weakly coupled oscillators exhibits synchronization in a parameter range. The synchronization has an anti-phase nature for an identical pair. However, it tends to be more in-phase for a non-identical pair with a rather weak coupling. A system of three identical oscillators (1, 2, and 3) coupled in a line (with two springs k12=k23) exhibits synchronization with two of them (1 and 2 or 2 and 3) being nearly in-phase. These collective behaviours are systematically estimated using the phase reduction method.


Introduction
Synchronization is ubiquitous in nature as there are numerous natural networks of nonlinear dynamical systems (Pikovsky et al., 2003). Because faults that cause earthquakes or seismogenic processes can be described as nonlinear dynamical systems, synchronization may occur in fault behavior (Scholz, 2010). The standard picture for the occurrence of interplate earthquakes is that a fault segment elastically driven by one plate, under the frictional resistance by another plate, exhibits a stick-slip motion that causes near-periodic spikes. A group of such segments can collectively cause recurring earthquakes with some statistical regularity (e.g., Scholz, 2002;Kawamura et al., 2012). Although many factors about the interaction between fault segments are still unknown, some evidence suggests that they can exhibit synchronization (Rubeis et al., 2010). For example, Chelidze et al. (2005) reported that a stick-slip object in a laboratory setting was entrained by a periodic force. Scholz Correspondence to: N. Sugiura (nsugiura@jamstec.go.jp) (2010) statistically determined that the occurrence of earthquakes in some regions was clustered. He reported that synchronous clusters of ruptures of several faults were identified in the south Iceland seismic zone, the central Nevada seismic belt, and the eastern California shear zone. Meanwhile, Mitsui and Hirahara (2004) successfully demonstrated that the numerically modeled coupled stick-slip oscillators exhibited some degree of synchronization. They used a simple spring-slider system composed of several mutually coupled stick-slip oscillators to capture the nature of the earthquake generation cycle along the Nankai trough, which is located in a zone of high seismicity where multiple segments that constitute the fault zone have been reported to rupture almost simultaneously (Ishibashi, 2004a). It is worth noting that they found that a pair of coupled oscillators with slightly different parameter sets synchronized even for weak coupling (Fig. 6 of Mitsui and Hirahara (2004)), although their emphasis was on cases with strong coupling between oscillators.
In spite of these observations, there has been little research that provides a specific description of the conditions for synchronization and how phases behave collectively. In this regard, we focus on the time evolution of the phases to elucidate the synchronization dynamics behind such collective behaviors and how phases are locked in the synchronization.
The occurrences of some earthquakes are nearly periodic (e.g., Matsuzawa et al., 2002;Ishibashi, 2004b;Sykes and Menke, 2006); thus, the generation process can be well modeled as a limit-cycle oscillation. The timing of a limit-cycle oscillation can be described by a single phase variable. If the limit-cycles are somehow connected, they should interact with each other and exhibit some collective behavior as a consequence of the attraction or repulsion between them in terms of the phase. The phase reduction method (Kuramoto, 1984) enables us to quantify the rate at which the progress of an oscillator phase is affected by another oscillator, thereby offering a powerful analytical tool to approximate the limit-cycle dynamics as a closed equation 2 SUGIURA ET AL.: SYNCHRONIZED OSCILLATORS for only a single phase variable.
We shall confine our attention to simple systems of only a few oscillators that remain close to a common limitcycle orbit, rather than the complicated ones that may produce chaotic motion (e.g., Turcotte, 1990, 1992;Abe and Kato, 2012), so that we can extract some regularity from the collective behavior of the oscillator system. This setting, of assuming almost homogeneous system of limit cycle oscillators, looks reasonable in the light of observations. In fact, there are some seismic zones that consist of fault segments that have quite similar recurrence periods. The deviation of the earthquake generation periods between different segments along the Nankai trough is a few years, much smaller than the periods themselves, ∼ 1 × 10 2 year (Ishibashi, 2004a). Likewise, Scholz (2010) points out that synchronization occurs within systems of evenly spaced, subparallel faults with very similar slip rates.
In this study, we quantitatively analyze how a single slider oscillates under the rate-and state-dependent friction against a plate motion using a bifurcation analysis and centermanifold reduction method. Then, we identify when and how coupled sliders driven by a plate synchronize as a collective substance using a phase reduction method.

The spring-slider-dashpot system
It is well established that a fault segment that can cause earthquakes is well described by a spring-slider system (e.g., Perfettini and Avouac, 2004) subjected to a rate-and statedependent friction (Dieterich, 1979;Ruina, 1983;Scholz, 1998); this model exhibits a limit cycle oscillation.
Our research interest, therefore, is in spring-coupled sliders ( Fig. 1) that are driven by a common plate through spring and dashpot arrangements set for each slider (e.g., Rice, 1993;Cochard and Madariaga, 1994), against the frictional resistance by another plate. The equations of motion for the i-th slider are where x i is the position of the slider, x 0 i and x 0 ij are the lengths of springs at rest, k i is the spring constant between the slider and plate, k ij is the spring constant between a pair of sliders, V i is the velocity, G is the rigidity, c is the shear wave velocity, and V p is the constant velocity of the plate. The frictional force τ i has a rate-and statedependent form that can be represented as (Ruina, 1983;Dieterich and Kilgore, 1994): where a i and b i are frictional parameters, σ i is the normal stress, V * and θ * are the arbitrary reference velocity and state, respectively, and µ * i is a reference frictional coefficient. The state variable θ i obeys an aging law represented by (Ruina, 1983;Linker and Dieterich, 1992): where L i is the characteristic length. Under a quasistatic approximation where the inertia m i d 2 x i /dt 2 is sufficiently small (Gu et al., 1984;Perfettini and Avouac, 2004;Perfettini et al., 2005;Kano et al., 2010Kano et al., , 2013, the governing equation for V i can be derived as where In accordance with the typical applications of the model to the seismogenic process, we assume that the parameters are in the range of We also assume that all of the initial states are placed in the first quadrant: In Section 3, we investigate the basic properties of a single oscillator. After introducing the phase reduction method in Section 4, we analyze the properties of synchronization, which occurs in a two-oscillator system, in Section 5. We mention some extensions to a three-oscillator system in Section 6.

The Dieterich-Ruina oscillator
Here, we investigate the basic properties of a single oscillator using the bifurcation and perturbation analyses.

Governing equations
Dropping the index i in Eqs. (4) and (5) for simplicity, we obtain the following equations describing a single oscillator: This is a two-dimensional dynamical system with six parameters (k,V p ,g,A,B,L). Hereafter, the dynamical system described by Eqs. (10) and (11) is called the Dieterich-Ruina SUGIURA ET AL.: SYNCHRONIZED OSCILLATORS 3 oscillator, and the state vector is denoted as X = (θ,V ) T . For the simplicity of analytical expressions, we use a new parameter set (µ,V p ,C,d,q,L) that is defined as where µ serves as a bifurcation parameter. Here, we investigate how the system behaves as µ changes. This system has a unique equilibrium point at X 0 = (L/V p ,V p ) T , which is given by the intersection of the nullclines: One of the important facts concerning the linear structure of the system around an equilibrium point is that the Jacobi matrix J has a characteristic sign pattern given by which represents a substrate-depletion system (Arcuri and Murray, 1986). The two eigenvalues of the Jacobi matrix at X 0 are

Stable spiral: µ < 0
The equilibrium point is a stable spiral when −2q/C < µ < 0 because the eigenvalues of J are a complex conjugate pair, and have a common negative real part.

Hopf bifurcation: µ = 0
At the very instance when µ = 0, the equilibrium point begins to lose its stability. The system encounters a Hopf bifurcation because the Jacobi matrix has a pair of imaginary eigenvalues The corresponding eigenvectors are and its complex conjugate, U . U and U span the plane containing linear solutions. By introducing a complex amplitude, W (t), the neutral solution of the system is expressed as where c.c. represents the complex conjugate. The graph containing the solution indicates an elliptic orbital motion, while the complex amplitude is an arbitrary complex constant at this stage (µ = 0) if we neglect nonlinear terms.

Weakly nonlinear: µ 0
When the bifurcation parameter µ becomes slightly larger than 0, the equilibrium point becomes an unstable spiral because the eigenvalues of J are a complex conjugate pair with a common positive real part. Here, we develop an analytical expression for the asymptotic solutions in a weakly nonlinear regime, by expanding Eqs. (10) and (11) as a Taylor series in terms of the deviation u ≡ X − X 0 (See appendix A), and using U in Eq. (20) and its dual, U * (a left eigenvector), which is given by: With the expansion and eigenvectors, we can compute the coefficients for a small-amplitude equation near the Hopf bifurcation following the center-manifold reduction method described in Kuramoto (1984). Assuming the solutions are in the form of Eq. (21), the time evolution of the complex amplitude can be described by the Stuart-Landau equation as This system encounters a supercritical bifurcation to a stable limit cycle, because the supercriticality condition Reβ > 0 is derived from 0 < d = gV p /(A + gV p ) < 1. Note that the type of the bifurcation may have some dependence on the laws of friction and assumptions made on the equation of motion (Gu et al., 1984;Putelat et al., 2010). In the original vector form, the limit-cycle solution of Eq. (24) is given by which graphically describes an elliptic orbital motion. The modulus, R s , and frequency shift,ω, are scaled with µ 1/2 and µ, respectively.

SUGIURA ET AL.: SYNCHRONIZED OSCILLATORS
We performed numerical integrations of Eqs. (10) and (11) to simulate the limit-cycle oscillation near the Hopf bifurcation point for three cases with µ = 10 −5 , 10 −4 , and 10 −3 Nm −2 . The time integrations were performed with the fourth-order Runge-Kutta scheme containing variable time step-sizes (Press et al., 1992). The rest of the parameters were set according to a previous study by Kano et al. (2010) for an inter-plate earthquake occurred on September 25, 2003 in Hokkaido, Japan: (V p ,g,A,B,L) = (3.17 × 10 −9 ms −1 , 5.00 × 10 6 Nm −3 s, 1.50 × 10 5 Nm −2 , 2.20 × 10 5 Nm −2 , 1.00 × 10 −2 m); these values also serve as the standard set of parameter values for this study. In Fig. 2, we show the results with the orbits of the limit cycle compared to those derived using Eq. (27). The corresponding orbits are in good agreement when µ is small.
In the context of seismogenic processes, the analytical solution (Eq. (27)) in the weakly nonlinear regime may offer a simplified description of slow earthquakes (e.g., Yoshida and Kato, 2003;Helmstetter and Shaw, 2009), which can be viewed as sustaining aseismic oscillations in which the slip instability is sufficiently weak (Kawamura et al., 2012). In particular, the frequencies in Eqs. (22) and (29) can be used to evaluate the recurrence intervals of such earthquakes.

Limit cycle: µ > 0
When we increase µ, the system will enter a strongly nonlinear regime. The equilibrium point becomes either an unstable spiral (when 0 < µ < 2q/C) or unstable node (when µ > 2q/C). Then, the Poincaré-Bendixson theorem (e.g., Strogatz, 2001) ensures the existence of a limit cycle within some region surrounding the equilibrium point, because we now have an unstable equilibrium point with a surrounding trapping region, R. Appendix B describes how flows are trapped into the region. Figure 3 shows an example of a limit cycle orbit derived by numerically integrating Eqs. (10) and (11). The orbit appears more polygonal than elliptical and extends over a wide range in the first quadrant.

The phase reduction method
Here, we introduce the phase reduction method for general limit-cycle oscillators, as well as its specific representation for weakly nonlinear oscillators.

Limit-cycle oscillators
A system of coupled self-sustained oscillators can be described by where we assume that the system dX/dt = F (X) behaves by itself as a limit-cycle oscillator and that the system de-scribed by Eq. (30) has an oscillatory behavior similar to it, including the frequency and orbit. Provided that the oscillators have similar properties and are weakly coupled, the phase reduction method (Kuramoto, 1984), shown below, is applicable to the system. Using the period, T , and the frequency, ω, for the limit cycle of the system dX/dt = F (X), we can define the phase, φ, of a state that is determined up to an integral multiple of T , which varies from 0 to 2π. The time evolution of the phase obeys where φ i is the phase of the oscillator i, δω i is the frequency deviation of oscillator i from the original limit cycle frequency, and Γ ij is the phase coupling function (hereafter, the PCF) between the oscillators i and j, which is periodic with a period of 2π. These terms are defined as the averaged values of the deviation terms in Eq. (30) over a period of the limit cycle under the action of phase sensitivity, Z(φ), (a row vector): Here, Z(φ) coincides with a left Floquet eigenvector, with eigenvalue 0, for the linearized equation around the limit cycle. Refer to Kuramoto (1984) for the details of the phase reduction method discussed here. This procedure is applicable to the system containing Dieterich-Ruina oscillators (Eqs. (4) and (5)), provided that both the parameter differences and coupling intensities of the oscillators are small enough to be treated as a perturbation. Substituting the specific functions in Eq. (5) into Eq. (33), we obtain the phase description of the system: where V * is the phase sensitivity for V . Note that V and V * are defined along a stable orbit of a single oscillator without coupling, which has a frequency ω.

Weakly nonlinear oscillators
Suppose we have a system of weakly nonlinear oscillators that are identical and mutually coupled. Near the Hopf bifurcation point, each oscillator can be described by Eq. (24) and a coupling term, which is supposed to be small: Normalizing the equations to t ′ = (µReα)t and By treating each oscillator as a two-dimensional system with independent variables (ReW ′ i ,ImW ′ i ) T , we can analytically derive the PCF for this complex Ginzburg-Landau-type equation (Kuramoto, 1984): For the case of the weakly nonlinear Dieterich-Ruina oscillators, (Eqs. (4) and (5) near the Hopf bifurcation point), the coupling coefficient, γ, is defined in the same manner as α in Eq. (25): Substituting Eqs. (25), (26), and (40) into Eq. (38), we get the coefficients in Eq. (37): Thus, the PCF (Eq. (39)) for the weakly nonlinear Dieterich-Ruina oscillator is characterized by In particular, the inequality (42) indicates that the coupling has an anti-phase nature (dΓ/dψ(0) > 0, dΓ/dψ(π) < 0)). Figure 4 shows the PCF as a function of the phase with the same parameters as in Fig. 2.

Two-oscillator system
Here, we explore when and how synchronization occurs in the system of two mutually coupled Dieterich-Ruina oscillators. We assume the two oscillators are identical except for a slight difference in the value of B i . To confirm the applicability of the phase reduction method to the stick-slip oscillator system, we examine the properties of the synchronization in two different ways. First, we observe the synchronization through numerical integrations of a coupled oscillator system. Second, we derive the PCF for the phase equations using the results from the numerical integration of a single oscillator system and its adjoint. Then, we determine some quantities from the plot.

Numerical Integrations
We performed numerical integrations of a discrete-time version of Eqs. (4) and (5), for a pair of coupled oscillators: a reference oscillator (oscillator 1) and a second oscillator (oscillator 2). Oscillator 1 had the following parameters: (k,V p ,g,A,B,L) = (1.00 × 10 5 Nm −3 , 3.17 × 10 −9 ms −1 , 5.00 × 10 6 Nm −3 s, 1.50 × 10 5 Nm −2 , 2.20 × 10 5 Nm −2 , 1.00 × 10 −2 m), and the natural frequency was 1.0687876 × 10 −9 s −1 . Two identical oscillators are coupled for case 0. For cases 1, 2, and 3, we used oscillator 2 that has the same set of parameters as oscillator 1 except for B = 2.2025 × 10 5 , 2.225 × 10 5 , and 2.25 × 10 5 Nm −2 , respectively. The natural frequencies of oscillator 2 in cases 1, 2, and 3 were 1.0652082 × 10 −9 s −1 , 1.0340304 × 10 −9 s −1 , and 1.00143857 × 10 −9 s −1 , respectively. We used a common coupling strength of K = k 12 = k 21 = 3×10 3 Nm −3 for all cases, based on the one used in Kano et al. (2010), which was derived through the inversion of strain rate from the GPS observation. We also checked that the values of B and K were within the range of application of the phase reduction method (See appendix C). The time integrations are performed using the same method described in Section 3.4. Figure 5 shows the results of case 0, in which the oscillators synchronize at the phase difference ψ = −3.14 (anti-phase). Figure 6 shows the results of case 1, in which the oscillators synchronize at the phase difference ψ = −1.18 (out-of-phase). Figure 7 shows the results of case 2, in which the oscillators synchronize at the phase difference ψ = −7.07 × 10 −3 (almost in-phase). Figure 8 shows the results of case 3, in which they exhibit no synchronization. These phase differences and synchronized oscillator frequencies are listed in Table 1.

Application of the PCF
In this setting, the evolution of the phases can be described as where the PCF is defined in Eq. (35), and the difference between the natural frequencies is estimated to be

SUGIURA ET AL.: SYNCHRONIZED OSCILLATORS
Taking the difference between Eqs. (44) and (45), we obtain the time evolution of the phase difference, ψ = φ 1 − φ 2 : Using a primitive function on the right-hand side of Eq. (48), we find that the phase difference obeys a gradient dynamical system As t → ∞, the state approaches a stable point at the bottom of the potential U . The realization of synchronization is equivalent to the existence of a phase difference ψ sync that satisfies Taking the average of Eqs. (44) and (45), we obtain the time evolution of the phase average, ϕ = (ψ 1 + ψ 2 )/2: where ω = ω + (δω 1 + δω 2 )/2. When synchronization is achieved, the frequency is shifted to We calculated the phase sensitivity, V * , with a relaxation method (Ermentrout, 1996;Ermentrout and Terman, 2010), using a numerical integration of the adjoint model of the Dieterich-Ruina oscillator. The integration is also performed using a fourth-order Runge-Kutta scheme with variable timestep sizes (Press et al., 1992). Figure 9 shows the phase sensitivity, or the values of V * , during a time interval. The value of the sensitivity remains positive for most of the period except at the moment when a slip event occurs. After the slip event, the sensitivity starts to increase for a while, following which it gradually decreases. Using the calculated values of V , θ, and V * as functions of phase, we have also calculated the PCF for the oscillator according to Eq. (35). Figure 10(a) shows the PCF as a function of phase. Figure 10(b) shows the PCF on negative and positive half-planes of phase in a logarithmic scale. The PCF is classified as an anti-phase type as in Fig. 4, although the shape does not resemble a sine curve.
By checking the positional relation between the horizontal line,Γ = −∆ω/(2K), and antisymmetric part,Γ a , of the PCF curve in Fig. 10, we can determine whether Eq. (52) subject to inequality (53) has a solution, i.e., we can evaluate whether synchronization is achieved. In this setting, synchronization is expected in the range of −6.5 × 10 −15 < −∆ω/(2K) < 6.5 × 10 −15 kg −1 m 2 s. If there is an intersection between the horizontal line and antisymmetric part,Γ a , of the PCF, in addition toΓ a being a decreasing function of the phase at that point, then the synchronization is achieved with the difference of phase at which the intersection is located, as indicated in Fig. 10(a). The frequency of synchronized oscillators is also derived using the symmetric part,Γ s , of the PCF according to Eq. (56).

Comparison of the results of numerical integration and phase reduction
In table 1, important quantities representing the synchronization properties are summarized: the difference of the natural frequencies ∆ω, phase difference ψ sync , and frequency dϕ/dt| sync . Data in the parentheses are the estimated values for the synchronization properties of the oscillator pairs, which are derived from the intersection of the PCF and a horizontal line. The estimated values from the PCF are in reasonable agreement with the corresponding ones by numerical integration. This indicates that the synchronization properties of coupled oscillators can be quantitatively estimated using the phase reduction method when the coupling is sufficiently weak. The PCFΓ a has a pretty complicated "microstructure" near |ψ| ≃ 0, a flat hill-like structure 0 ψ < 10 −2 with a sudden jump to the origin, as shown in Fig. 10(b). Thus, we cannot decide the exact phase difference at which the oscillators are nearly synchronized in-phase. However, we are sure that the phases never become exactly in-phase, whereΓ ′ a violates the inequality (53). According to the phase reduction method, the range of parameters in which two oscillators synchronize is estimated to be where ∆ represents the difference between two oscillators. This gives |∆B| < 2.7 × 10 3 Nm −2 for the parameters used here, which is consistent with the results of numerical integrations.

Three-oscillator system
Here, we extend the analysis to a system of three identical mutually coupled Dieterich-Ruina oscillators. Each oscillator in the system is assumed to be described by Eqs. (4) and (5) and contain the same parameter set as oscillator 1 in Section 5. We consider two different coupling topologies: a periodical coupling in a ring with spring constants k 12 = k 23 = k 31 = K and a non-periodical coupling in a line with k 12 = k 23 = K, k 31 = 0. In terms of phase, the state of the three-oscillator system can be characterized by the phase differences between the oscillators, ψ 1 = φ 1 − φ 2 and ψ 3 = φ 3 − φ 2 .

Numerical Integrations
We performed numerical integrations of a discrete-time version of the original system of differential equations for the two types of coupling patterns. Figure 11 shows the results for the periodical coupling. After convergence, the three oscillators share a common phase difference of 2π/3, i.e., (ψ 1 ,ψ 3 ) ≃ ( 2 3 π,− 2 3 π). Figure 12 shows the results for the non-periodical coupling. Although the convergence is rather slow, the oscillators gradually synchronize at phase differences near (ψ 1 ,ψ 3 ) ≃ (7.0 × 10 −3 ,2.6).

Application of the PCF
By applying the phase reduction method to the system of three identical oscillators, the evolution of the phases can be described as whereΓ = Γ/K. Differences between the three equations in (58) give the time evolution for the phase differences as in Eq. (48). The time evolution of the system of periodically coupled oscillators can be written as The time evolution of the system of non-periodically coupled oscillators can be written as Each three-oscillator system is thereby reduced to a twodimensional dynamical system for the phase differences (e.g., Aihara et al., 2011); this two-dimensional system has a symmetry because ψ 1 and ψ 3 are interchangeable. Similar to conditions (52) and (53) for a system of two oscillators, the synchronization of the three-oscillator system is expected to be realized at the stable equilibrium points of the phase flow in the (ψ 1 ,ψ 3 )-plane; these equilibrium points emerge as intersections of the nullclines for the phase flows (Figs. 13  and 14). In the upper-right part of Fig. 14, the nullclines for dψ 1 /dt = 0 and dψ 3 /dt = 0 nearly overlap because Γ(ψ 1 ) is almost equal to Γ(ψ 3 ) owing to a rather flat region of Γ ( Fig. 10(b)) they share in this range.

Comparison of the results of numerical integration and phase reduction
The triphase synchronization (e.g., Aihara et al., 2011) in the periodically coupled system (Fig. 11) is achieved because the phase oscillators exclude each other with an equal intensity owing to the anti-phase nature of the PCF (Fig. 10). It corresponds to a stable spiral in the fourth quadrant of the phase plane (Fig. 13). The synchronization in the non-periodically coupled system (Fig. 12) corresponds to one of the two stable nodes in the first quadrant of the phase plane (Fig. 14). The reason for the slow convergence for the latter case is that the orbit of the phase differences should follow a static pathway along one of nearly overlapped nullclines mentioned above.
In each three-oscillator system, the phase flow has a pair of stable equilibrium points at a symmetric position in the phase plane with different basins of attraction. Hence, the convergence of the phase differences is dependent on which basin the initial condition belongs.

Conclusions
The Dieterich-Ruina oscillator can be viewed as a selfsustained oscillatory system with two degrees of freedom. This concisely describes the stick-slip motion of a slider driven by a plate through a spring and dashpot against a rateand state-dependent friction. When the bifurcation parameter µ = σ(b − a) − GV p /(2c) − Lk passes through zero, it encounters a supercritical Hopf bifurcation, and an asymptotic analytical solution (Eqs. (22), (27)-(29)) in the weakly nonlinear regime is available for µ 0, which may serve as a formula for evaluating the recurrence intervals of slow earthquakes if the slip instability is sufficiently weak.
Some collective behaviors are found for a pair of weakly coupled Dieterich-Ruina oscillators. The two-oscillator system that has similar weakly coupled oscillators exhibited synchronization for some combinations of the coupling strength and similarity of the oscillators. Synchronization is expected in the parameter range of inequality (57). Even though different systems of oscillators should have different criteria, a simple model for the earthquake generation cycle along the Nankai trough exhibited synchronization in a similar range of η < 0.35 (Cases 1, 2, and 3 of Mitsui and Hirahara (2004)), which suggests that synchronization can occur in seismogenic process.
The synchronization is anti-phase for an identical pair; however, their phases tend to align for non-identical pairs with weak coupling. The phase behavior was quantitatively estimated using the phase coupling function for the oscilla-8 SUGIURA ET AL.: SYNCHRONIZED OSCILLATORS tor. It is interesting that a pair of non-identical oscillators with weak coupling can nearly cause an in-phase synchronization. This suggests the possibility of sequential occurrences of adjacent earthquakes.
Distinct phase alignment behaviors were found for threeoscillator systems. The system of three identical oscillators equally coupled in a ring exhibits a triphase synchronization, in which they arrange themselves such that they are out-ofphase with respect to each other by 2π/3. In contrast, if three identical oscillators (1, 2, and 3) are equally coupled in a line with spring constants k 12 = k 23 , k 31 = 0, then oscillators 1 and 2 or oscillators 2 and 3 become nearly in-phase, while the other remains nearly anti-phase. The synchronization properties were quantitatively estimated using the phase reduction method.
These results demonstrate that synchronization should occur between several coupled oscillators in stick-slip motion, for which we can systematically use the phase reduction method as an analytical tool. In the context of seismogenic processes, the phase reduction method can be applied to the analysis of observed synchronization in a seismogenic zone that is presumed to consist of neighboring groups of faults moving at similar slip rates with mutual stress coupling (Scholz, 2010). Moreover, the method is still applicable even if the inertia is included in the model or another friction constitutive law is adopted. It may be of interest to examine how the phase coupling function changes its property according to these details of the modeling. In particular, it will be meaningful to specify the extent to which the inertia term will affect the timing of slip events.
The phase description (34) has a general form applicable to a system with an arbitrarily large number of the oscillators described by Eqs. (4) and (5), as long as the oscillators are weakly coupled. As a consequence of the anti-phase nature of the oscillator, which is evident from the inequality (42) or from the shape of the PCF (Fig. 4 or 10), an irregular pattern may emerge even in a homogeneous system with a large population of diffusively coupled oscillators. This is where we will be able to find the Benjamin-Feir instability developing phase turbulence (Kuramoto, 1984).

Appendix A The Taylor expansion
The Taylor expansion of Eqs. (10) and (11) in terms of the deviation u ≡ X − X 0 is as follows.
where h.o.t . denotes higher order terms.

Appendix B The trapping region
The region R can be constructed by bounding it with a hexagon H = ABCDEF in a logθ-logV plane, where A = (logθ 1 ,logV 1 ), B = (logθ 4 ,logV 1 ), C = (logθ 4 ,logV 3 ), D = (logθ 2 ,logV 2 ), E = (logθ 3 ,logV 2 ), and F = (logθ 5 ,logV p ). Using small positive values ǫ i , 1 ≤ i ≤ 6, we can define the constants for these positional coordinates as SUGIURA ET AL.: SYNCHRONIZED OSCILLATORS 9 An example of the trapping region is illustrated in Fig. 15. If we assign appropriate values to ǫ i , then all the trajectories in R will be confined within it. To be specific, we can set the diagonal segments CD, EF , and F A to be sufficiently steep, or vertical, such that any flows on them would be trapped. This is derived as follows. Slopes of the flows on a logθ-logV plane are defined as Here, we assess this quantity especially on the diagonal segments CD, EF , and F A.
-Segment CD Since CD has no intersections with nullcline I, and is placed on the upper right side of it, the quantity 1/(1 − V θ/L) − 1 has a finite negative value. Hence, if we assign V a large value, |γ| can be arbitrarily small. In other words, if we place the segment CD in a large V region, then the flows on it should have sufficiently gentle, or horizontal, slopes to be trapped in the region R.
-Segment EF On this segment, using V p ≤ V ≤ V 2 and 0 < V θ/L ≤ ǫ 5 < 1, we find From these three inequalities, we get an estimation for the negative slope γ: The rightmost inequality is a consequence of µ > 0. If we use ξ 1 as the slope of the segment EF , then flows on it should have sufficiently gentle slopes to be trapped in the region R.
-Segment F A On this segment, using V 1 ≤ V ≤ V p and 0 < V θ/L < 1, we find From these three inequalities, we get an estimation: Using the same method as the EF case, if we use ξ 2 as the slope of the segment F A, then flows on it should have gentle slopes to be trapped in the region R.
See Strogatz (2001) for the construction of trapping regions.

Appendix C The range of application of the phase reduction method
To investigate the range of application of the phase reduction method, we quantify here the weakness of heterogeneity and interaction of the oscillators in terms of the system of mutually coupled Dieterich-Ruina oscillators. Since the orbit of oscillator is well captured on a logarithmic scale as in Fig. 3, it is convenient to deal with the logarithm of the variables in this discussion. The time evolution of (logθ i ,logṼ i ) can be written in a dimensionless form: We can apply the phase reduction method if the perturbations caused by the oscillator difference and the coupling term are sufficiently smaller than the absolute value of the Floquet exponent for the amplitude mode of the limit cycle dX/dt = F (X). This condition ensures that the orbits of coupled oscillators stay in the neighborhood of the original limit cycle orbit owing to the restoring effect. The Floquet exponent for the amplitude mode of oscillator 1 in section 5 is estimated to be λ = −8.8 × 10 −10 s −1 , while the averaged perturbations caused by the oscillator difference and the coupling term are estimated to be where ∆B = B i − B j , K = k ij , and the integrations are performed along the limit cycle orbit. Substituting these into ∆λ h ,∆λ c ≪ |λ|, we get the conditions for ∆B and K: Furthermore, we can apply averaging over a period to derive the phase shifts if the averaged perturbation on the phase does not alter the natural frequency substantially. The natural frequency of oscillator 1 in section 5 is ω = 1.0 × 10 −9 s −1 , while the perturbations on the phase caused by the oscillator difference and the coupling term are estimated, using Eqs. (35) and (47), to be Substituting these into ∆ω h ,∆ω c ≪ ω, we get the conditions for ∆B and K: Taking into account these criteria, we choose these parameters in the range of 0 ≤ ∆B ≤ 5 × 10 3 Nm −2 and 0 ≤ K ≤ 3 × 10 3 Nm −3 . We have also checked directly that each orbit in the numerical integrations stays in the neighborhood of the original limit cycle orbit. Figure 16 shows a comparison of the orbits.  (Burridge and Knopoff, 1967), except that it is also equipped with dashpots and the friction on the bottom of the sliders is rateand state-dependent.  (1,1). The values of bifurcation parameter µ are set to 1 × 10 −5 Nm −2 (solid curves), 1 × 10 −4 Nm −2 (dashed curves), and 1 × 10 −3 Nm −2 (dotted curves). The values of k used here are obtained by setting k = (B − A − gVp − µ)/L and using the corresponding values of µ. The rest of parameters are set to (Vp,g,A,B,L) = (3.17 × 10 −9 ms −1 , 5.00 × 10 6 Nm −3 s, 1.50 × 10 5 Nm −2 , 2.20 × 10 5 Nm −2 , 1.00 × 10 −2 m).  Phase sensitivity, V * , as a function of time in years (red curve) calculated numerically using the relaxation method. The parameters are set to (k,Vp,g,A,B,L) = (1.00 × 10 5 Nm −3 , 3.17 × 10 −9 ms −1 , 5.00 × 10 6 Nm −3 s, 1.50 × 10 5 Nm −2 , 2.20 × 10 5 Nm −2 , 1.00 × 10 −2 m). The green curve is for V /Vp in a logarithmic scale. Note that V * becomes negative for some time periods in which V /Vp is large. The blue, green, and red curves are the antisymmetric partΓa, symmetric partΓs, and totalΓ, respectively. For each case in Table 1, the phase difference ψ of synchronized oscillators and the corresponding value of −∆ω/(2K) are indicated by a filled circle and an arrow, respectively.  11. The time evolution of V /Vp for three identical oscillators with a periodic coupling in a logarithmic scale. Red, blue, and green curves correspond to oscillator 1, 2, and 3, respectively. The variation from 5000 to 65000yr is not shown. The parameters are B = 2.20×10 5 Nm −2 for all oscillators, and k12 = k23 = k31 = 3 × 10 3 Nm −3 (periodic one-dimensional coupling). The three oscillators synchronize at the phase differences (ψ1,ψ3) ≃ ( 2 3 π,− 2 3 π). This synchronization corresponds to a stable spiral in Fig. 13. Red, blue, and green curves correspond to oscillator 1, 2, and 3, respectively. The variation from 5000 to 13825000yr is not shown. The initial states for three oscillators are different from each other. The parameters are B = 2.20 × 10 5 Nm −2 for all oscillators, and k12 = k23 = 3×10 3 Nm −3 ,k31 = 0 (non-periodic one-dimensional coupling). They synchronize at the phase differences of (ψ1,ψ3) ≃ (7.0×10 −3 ,2.6), where oscillators 1 and 2 are nearly in-phase. This synchronization corresponds to a stable node in Fig. 14 Fig. 13. Flow directions and nullclines on the (ψ1,ψ3)-plane for the phase flow of three identical oscillators that are periodically coupled. Arrows indicate the flow direction. Green and red curves represent the nullcline dψ1/dt = 0 and dψ3/dt = 0, respectively. Stable spirals are located at ± 2 3 π,∓ 2 3 π . The origin is an unstable node, and the three saddles are around ( Flow direction and nullclines on the (ψ1,ψ3)-plane in a logarithmic scale for the phase flow of three identical oscillators that are non-periodically coupled. Arrows indicate the flow direction. Green and red curves represent the nullcline dψ1/dt = 0 and dψ3/dt = 0, respectively. Stable nodes are located around 8 × 10 −3 ,2.6 and 2.6,8 × 10 −3 . The origin is an unstable node, and saddles are around   . 16. A comparison of the orbits, graphs of (θVp/L,V /Vp), in a double logarithmic plane. The red curve, green dots, and blue triangles are the orbits of the original limit cycle, oscillator 1 and 2 in case 3, respectively. 3.62 × 10 −11 3.44 × 10 −11 1.76 × 10 −11 -(3.61 × 10 −11 ) (3.44 × 10 −11 ) (1.81 × 10 −11 ) ( -)