Earthquake Sequencing : Chimera States with Kuramoto Model 2 Dynamics on Directed Graphs

Recommendation to the Editor Does the paper contain new and significant results? Yes No Is the paper of an international standard? Yes No Is the presentation clear and concise? Yes No Does the paper put the obtained results into context, with relevant references? Yes No Is the length of the paper appropriate? Yes No Is the text fluent and precise? Yes No Are the title and the abstract pertinent and understandable to a wide audience? Yes No


Introduction
Earthquakes of differing magnitudes occur at different locations and depths in many tectonically active regions of the earth.The magnitude is the most widely used and theoretically studied earthquake parameter (Kanamori and Anderson, 1975;Hanks and Kanamori, 1979).The moment magnitude scale, M W , provides an estimate for all medium to large earthquake magnitudes.Continuous recording and analysis of earthquakes that occur in different regions of the earth have led to earthquake catalogues.These catalogues carry information about the epicenter and the estimated hypocenter, the time and the magnitude of the earthquakes, leading to a set of empirical rules for different earthquake regions and the global seismicity (Omori, 1895;Gutenberg and Richter, 1954;Bath, 1965;Bufe and Varnes, 1993;Utsu et al., 1995;Ogata, 2011).The empirical rules allow us to understand and expand on the inter-relationships between the earthquake magnitude and the frequency of occurrence of events, and the main shocks and their aftershocks in space and in time.
The earthquake catalogues have recently become the basis for Markov chain models of earthquake sequencing to explore probabilistic forecasting from the point of view of seismic hazard analysis (Nava et al., 2005;Cavers and Vasudevan, 2015).Cavers and Vasudevan (2015) have incorporated the spatio-temporal complexity of the earthquake recurrences (Davidsen et al., 2008;Vasudevan et al., 2010) into their Markov chain model.Intrinsic to earthquake sequencing studies is the observation made on scaling behavior and earthquake cycles (Turcotte, 1997;Rundle et al., 2002Rundle et al., , 2003)).In this regard, fractal and fractal-rate stochastic point processes were found to be useful (Thurner et al., 1997).Telesca et al. (2011) applied such models to earthquake sequencing.Vasudevan and Cavers (2013) have recently extended the application of this model to study time-correlative behavior in earthquake sequencing by carrying out Fano factor and Allan factor analysis of a time series of state-to-state transition frequencies of a Markov chain.
One aspect of earthquake sequencing that requires a close look is a model for the non-linear dynamics of earthquakes.In this paper, we investigate the synchronization behavior of weakly coupled "earthquake oscillations".Such oscillations in the earth's crust and the epileptic brain show certain commonalities in that the distributions of energies and recurrence times exhibit similar power-law behavior (Herz and Hopfield, 1995;Rundle et al., 2003;Osorio et al., 2010;Chialvo, 2010).A growing interest in understanding the behavior of earthquakes and epileptic seizures with a view to exploring possible forecasting methods is one reason for the present study.In the case of epileptic seizures, the non-linear dynamics of pulse-coupled neuronal oscillations as an alternative to the Kuramoto (1975) model are under close scrutiny (Rothkegel and Lehnertz, 2014).To our knowledge, neither a simple Kuramoto model nor a modification of it has been worked out for earthquake sequencing studies.Mirollo and Strogatz (1990), Kuramoto (1991) and Rothkegel and Lehnertz (2014) considered the synchronization of pulse-coupled oscillators in which single oscillators release energy rapidly when they reach a trigger threshold and become quiescent for some time until they reach the trigger threshold again.Examples falling into this category are earthquakes and spiking neuronal activities (Herz and Hopfield, 1995;Beggs and Plenz, 2003;Rundle et al., 2002Rundle et al., , 2003;;Scholz, 2010;Karsai et al., 2012;Rothkegel and Lehnertz, 2014).Herz and Hopfield (1995) studied the collective oscillations with pulsecoupled threshold elements on a fault system to capture the earthquake processes.There are two timescales: the first is given by the fault dynamics defining the duration of the earthquake, and the second timescale is given by the recurrence time between "characteristic events", the largest earthquakes on a fault.The known recurrence times on several fault systems are 6 to 8 orders of magnitude longer than the duration of single events.Rundle et al. (2002) examined the selforganization in "leaky" threshold systems such as networks of earthquake faults.In their paper, they argued that on the "macroscopic" scale of regional earthquake fault systems, self-organization leads to the appearance of phase dynamics and a state vector whose rotations would characterize the evolution of earthquake activity in the system.Scholz (2010) invoked the Kuramoto model to represent the fault interactions, although no numerical synchronization-simulation results were presented.He postulated that the common oc-currence of triggering of a large earthquake by other earthquakes on nearby faults and the observation of space-time clustering of large earthquakes in the paleoseismic record were both indicators of synchronization occurring between faults.However, we need to bear in mind here that incorporating fault-fault interactions on a global scale involving all the networks of earthquake faults is formidable and nearly impossible.In this paper, we modify the simple non-linear mathematical model, the Kuramoto model with a phase lag, for the sequencing of global earthquake data.We show here that the solutions to the Kuramoto model with phase lag and with non-local coupling effects reveal the co-existence of synchronized and asynchronized states or chimera states for certain parameter values.We use this model as a precursor to our planned studies on other mathematical models such as integrate and fire models.
As alluded to earlier, there is a quiescence period between earthquakes in an earthquake zone, also known as the recurrence times.Since the globally recorded earthquake data are only available for a short time period, incorporating the recurrence times into the earthquake catalogue is impossible.Here, we consider the model proposed by Davidsen et al. (2008) to include the spatio-temporal complexity of recurrences by identifying the earthquakes occurring in close proximity to any occurred event in the record-breaking sequence.In this paper, we also investigate the Kuramoto model with a phase lag for the sequencing of global earthquake data influenced by the recurrences to point out the emergence of chimera states under certain conditions.

Mathematical model of the earthquake sequencing
The Kuramoto (1975) model for a large number of weakly coupled oscillators has become a standard template in nonlinear dynamical studies, pertinent to synchronization behavior, following the ground-breaking study of Winfree (1967).To apply this model to earthquake sequencing studies, we need to make a few justifiable assumptions and to incorporate certain essential features of earthquakes that we have come to know.For example, plate motions and, hence, plate tectonics (Stein, 1993;Kagan et al., 2010;DeMets et al., 2010) suggest that most of the earthquakes occur in and around plate boundaries because of the varying plate motions of the plates that uniquely encompass the earth's crust.In particular, different plates move at different rates and along different orientations, resulting in stress-field changes at the plate boundaries.When stress-field accumulation reaches, at a particular location or in a zone, a certain critical threshold, energy is released in the form of an earthquake.The relaxed system goes through the stress-build-up process again, a similar mechanism being operative in neuronal communication dynamics.We assume that there is a uniform stress increase during the quiescent period.Collective synchronization of threshold-coupled or pulse-coupled oscillators would be a candidate for such a study (Mirollo and Strogatz, 1990;Kuramoto, 1991;Rothkegel and Lehnertz, 2014).However, we defer the extension of their approach to earthquake sequencing studies to a future date.Since the quiescence period is 6 to 8 orders of magnitude longer than the event time duration, it would be an ideal platform on which to carry out this study.We surmise that the behavior of earthquake cycles noted in earthquake sequencing does not lend support, however, to a full synchronization or full asynchronization as a solution to this non-dynamics problem.One proven modification is the inclusion of non-local effects of the geometry of the system that has been shown to lead to a co-existence of partially synchronized and partially asynchronized states of oscillators as a steady-state solution.Such states, addressed as chimera states, are the subject of recent theoretical and experimental studies (Kuramoto and Battogtokh, 2002;Abrams and Strogatz, 2004;Abrams et al., 2008;Ko and Ermentrout, 2008;Omel'chenko et al., 2008;Sethia et al., 2008;Sheeba et al., 2009;Laing, 2009a, b;Laing et al., 2012;Martens et al., 2013;Yao et al., 2013;Rothkegel and Lehnertz, 2014;Kapitaniak et al., 2014;Pazó and Montbrió, 2014;Panaggio and Abrams, 2014;Zhu et al., 2014;Gupta et al., 2014;Vasudevan and Cavers, 2014a, b).We focus our present study on defining a Kuramoto model with a phase lag that would accommodate the existence of chimera states.The Kuramoto model has been extensively studied for a system made up of a large number of weakly coupled oscillators, where most of the physical problems are finite and can be described as nonlinear dynamics on complex networks (Acebrón et al., 2005;Arenas et al., 2008).In the realm of graph theory, complex networks can be cast as either undirected or directed graphs.In our studies on earthquake sequencing, we consider a directed graph as a representation of an earthquake complex network.The occurrence of chimera states as solutions to non-linear dynamics on both undirected and directed graphs has recently been investigated (Zhu et al., 2014;Vasudevan and Cavers, 2014a).As a precursor to studying earthquake sequencing with real data from the earthquake catalogues, we investigated the Kuramoto model on synthetic networks that mimic Erdös-Rényi random networks, small-world networks, and scale-free networks and directed graphs adapted from them, and examined chimera-state solutions (Vasudevan and Cavers, 2014a).For the earthquake sequencing studies here, we use the following Kuramoto model with a phase lag, α, with non-local coupling effects terms added explicitly: Here, θi is the time derivative of the phase of the ith oscillator.The angle α (0 ≤ α ≤ π/2) corresponds to the phase lag between oscillators i and j .G ij is the non-local coupling function that depends on the shortest path length, d ij , between oscillators i and j in the complex network: (2) K is the global coupling strength and κ is the strength of the non-local coupling.For convenience, we use a constant natural frequency for all the oscillators, i.e. homogeneous case, and, thus, we could use ω i = 0 for i = 1, . . ., N.Although we have not investigated the influence of the global coupling strength on the steady-state solution of the Kuramoto model, we treat this term as constant, in particular K = 1, based on observations made by Zhu et al. (2014).
We would like to stress that the model in Eq. ( 1) is not a pulse-coupled or threshold-coupled oscillator model.Although it would be appropriate to consider a variation of the Kuramoto model such as the Shinomoto-Kuramoto model (Shinomoto and Kuramoto, 1986;Sakaguchi et al., 1988;Lindner et al., 2004), we limit ourselves to a simpler model that does not include the excitable behavior of the model.We intend to use this model as a precursor to our planned studies on other mathematical models such as integrate and fire models.
A comment on the phase-lag parameter, α, in Eq. ( 1) is also in order.Panaggio and Abrams (2014) interpret the phase lag as an approximation for a time-delayed coupling when the delay is small.The value of α used is (π/2) − 0.10.In some ways, we treat the phase lag as a proxy for time delay.As Panaggio and Abrams (2014) demonstrate in their paper, the value of α determines a balance between order and disorder.We have not done an exhaustive search on the α parameter for the cases discussed in this paper.
Here, we construct a directed graph of earthquake events from the Incorporated Institutions for Seismology (IRIS) earthquake catalogue for the time period between 1970 and 2014.We consider earthquake events with magnitudes exceeding or equal to M W = 5.5 observed to a depth of 70 km.We partition the general latitude-longitude map of the earthquake events into a grid.We show two maps of such grid matrices (Fig. 1).A cell in a smaller grid (128 × 128) could have higher multiplicity of earthquake events than that in a large grid (1024 × 1024).We consider the coordinates of the topological center of each cell to represent the coordinates of the earthquake events that fall into that cell.Thus, we explore the effect of hubs and community effects by looking at transition probability matrices generated from grids of different orders such as 128, 192, 256, 512, and 1024 representing the seismicity map on a global longitude-latitude grid (Table 1).We compute the transition probability matrix and the shortest-path distance matrix for the directed graphs resulting from the catalogue considered.To keep the Kuramoto model simple, we assume a constant phase lag, α, in the phase of the ensemble of oscillators.The value of α used is (π/2) − 0.10.We relax this condition in subsequent simulations.The most difficult parameter to deal with here is the period of quiescence after the energy release following a certain stress threshold.We incorporate the build-up of the threshold effect indirectly by positing the inclusion of earthquake recurrences  in transition probability matrices.Here, we use the spatiotemporal recurrences based on the record-breaking model of Davidsen et al. (2008).In all our initial simulations, we ignore the influence of amplitude effects on the stability of the chimera states.We carry out simulations on the Kuramoto model for 200 000 time steps for the 128 × 128 oscillator grid matrices and for the 1024 × 1024 oscillator grid matrices.We report here the preliminary results of our simulations.

Results
We report the Kuramoto model experimental results for oscillators resulting from 128 × 128, 192 × 192, 256 × 256, 512 × 512, and 1024 × 1024 grids of the latitude-longitude map of the earthquakes.We consider a total of 13 190 earthquakes.We construct the transition probability and the shortest-path distance matrices for the grids without ("nonrecurrence" results) and with the consideration of the spatiotemporally complex recurrences ("recurrence" results), as shown in Figs. 2 and 3.
To represent the results, we use snapshots of three attributes (Zhu et al., 2014): (i) the phase profile, (ii) the effective angular velocities of oscillators and (iii) the fluctua-tion of the instantaneous angular velocity of oscillators.The effective angular velocity of oscillator i is defined as (3) Here, we take T = 1000 so that the effective angular velocities of the oscillators are averaged over the last 1000 time steps.We take t 0 = 199 001 for the 128 × 128 grid and for the 1024 × 1024 grid.
The fluctuation of the instantaneous angular velocity, σ i , of an oscillator i around its effective velocity is defined as If σ i = 0, then oscillator i rotates at a constant angular velocity.We show the non-recurrence and recurrence results obtained from the behavior of the last 1000 time steps of the simulations involving 200 000 time steps.We present the results for the 128 × 128 grid without and with recurrences in Figs. 4 and 5 for the three attributes using κ = 0.10.Figures 6  and 7 show these attributes for the 1024 × 1024 grid without and with recurrences, respectively, for κ = 0.1.Whether or not the Kuramoto model reaches the steady state, we examine the ratio of the number of coherent or synchronous oscillators to the total number of oscillators or "chimera index" as a function of the number of time steps.Here, we carry out 200 000 time steps.After every 20 000 time steps, we look at the chimera index for the last 1000 time steps.As an example, in Fig. 8, we find the asymptotic behavior of the scatter of the chimera index for ten such intervals for the 128 × 128 grid for κ = 0.10, suggesting that the Kuramoto model has reached the steady state.We investigate the influence of the non-local coupling coefficient, κ, on the chimera index for each grid and summarize our results for the non-recurrent 128 × 128 and 1024 × 1024 grids in Fig. 9.
Most of the initial computations reported in this work were on a HP C7000 chassis cluster system with dual-core 2.4 GHz AMD Opetron processors at the high-performance computing facility at the University of Calgary.We carried out a series of runs for 200 000 time steps on a Mac Pro Six-Core Intel Xeon E5 3.5 GHz, 16 GB RAM desktop work station and on a Dell PowerEdge R910 with Intel Xeon E7-4870 2.40 GHz 256 GB RAM processors, and we used the Matlab ODE113 solver to solve the Kuramoto model.

Building the directed graph
Earthquake sequencing is a well-studied problem in earthquake seismological communities around the globe, and yet it hides a suite of phenomenological mysteries that stand in the way of successful earthquake forecasting.One of the first steps in carrying out any investigative work on earthquake sequencing is to look at the global seismicity map such as the one posted by IRIS on a regular basis, with continuous updating of the associated catalogue.In Fig. 1a and b, we summarize the cumulative results of the catalogue for magnitudes of earthquakes exceeding M w = 5.5 and the depths of occurrence not exceeding 70 km, recorded between January 1970 and September 2014.One difference in the two figures lies in the coarseness of the gridding, with the first one being coarser than the second.A cursory glance at the figures immediately suggests the relevance of plate tectonics in that most earthquakes seem to occur at and around plate boundaries.A broad classification of these earthquakes could consist of the following categories: strike-slip earthquakes, subduction-zone seismicity, oceanic earthquakes, continental extensional regimes, intraplate earthquakes, and slow earthquakes (Scholz, 2002).The interplay between these remains a topic of research among seismologists.In general, fault systems play an important role in understanding the cause and recurrence of earthquakes.Scholz (2002) provides an excellent account of the mechanics of earthquakes and faulting.Ben-Zion and Sammis (2003) examined the continuum-Euclidean, granular, and fractal views of the geometrical, mechanical, and mathematical nature of faults and concluded that many aspects of the observed spatio-temporal complexity of earthquakes and faults might be explained using the continuum-Euclidean model.They contended that a continuum-based description would provide a long-term attractor for structural evolution of fault zones at all scales.The underpinning point in these works is the importance of the faulting in earthquake processes.Earthquakes are known to occur at different depths.Excepting in instances where there are surface ruptures as a result of earthquakes, fault zones at seismogenic depths in kilometers cannot be directly observed (Ben-Zion and Sammis, 2003).Continued geological mapping and high-resolution geophysical measurements afford a mechanism to improve our understanding of the fault zones.Rundle et al. (2003) took a statistical physics approach in emphasizing the significance of faults and fault systems as high-dimensional non-linear dynamical systems characterized by a wide range of scales in both space and time, from centimeters to thousands of kilometers, and from seconds to many thousands of years.The signature of the residual behavior in these systems is chaotic and complex.Understanding the coupling between different space scales and timescales to comprehend the non-linear dynamics of the fault systems is not an easy problem.In this regard, any attempt to explore the possibilities that accrue from non-linear dynamics studies is welcome.
In earlier studies on model and theoretical seismicity (Burridge and Knopoff, 1967;Vieira, 1999), special attention was paid to finding out whether chaos was present in the symmetric non-linear two-block Burridge-Knopff model for earthquakes.Vieira (1999) demonstrated with a three-block system the appearance of synchronized chaos.A consequence of this study was the speculation that earthquake faults, which are generally coupled through the elastic media in the earth's crust, could in principle synchronize even when they have an irregular chaotic dynamics (Vieira, 1999).Going one step further would be to suggest that the occurrence of earthquakes and the space-timescale patterns they leave behind is a sound proxy for modeling and theoretical studies of the fault systems.It is this point that is pursued in this work.
In this study, we focus on the non-linear dynamics of weakly coupled oscillators.Each oscillator (corresponding to the occurrence of an earthquake) is a proxy for a fault system or network with known information on its location, the time when the earthquake event occurred, and magnitude.This defines an element in the earthquake sequence.A continued sequence of events is represented as a directed graph (Vasudevan et al., 2010;Cavers and Vasudevan, 2014;Vasudevan and Cavers, 2014b) with the vertices representing the earthquakes (and their attributes) and the arcs the connecting links between neighbors in a sequence.Figures 2a  and 3a show the transition matrices for the directed graphs of the two grids, 128 × 128 and 1024 × 1024 grids.The oscillator index is determined by the grid partition with nonzero cells labelled in row-by-row order.A log(log) display scale is used to highlight the "clustering".The level of clustering along the first leading off-diagonal elements of the transition matrix is highlighted and indicates the partition-  ing and the relative significance of the seismicity zones in the globe.However, this does not invoke any causality argument.Since the multiplicity of the earthquakes in the cells of the grids used varies from "zero" to a large number, for the reason mentioned concerning the Euclidean geometry mentioned earlier, inter-cell and intra-cell transitions populate the transition matrices.These transition matrices are not symmetric.The non-linear dynamics of weakly coupled oscillators on such matrices has not been fully understood.
As mentioned earlier, the quiescence period between earthquakes in an earthquake zone is what we interpret here as a recurrence period.Studies on plate-boundary motions (Bird, 2003;DeMets et al., 2010;Stein, 1993) will provide an insight into the recurrence period for earthquakes in certain major fault zones.Even in instances where knowledge of the recurrence periods is known, it is usually punctuated by random fluctuations, the statistics of which are not unknown.The quiescence period is analogous to the process in human brains that precedes epileptic seizures (Berg et al., 2006;Rothkegel and Lehnertz, 2014), the structure of which has been modeled using pulse-coupled phase oscillators.Such pulse coupling or threshold coupling remains to be quantified for earthquakes.We defer this aspect of the work to future studies.Furthermore, the historical seismicity data set is short and, therefore, any information to be drawn from global records will be insufficient.However, the recurrence model introduced by Davidsen et al. (2008) offers a simple remedy to the problem identifying the earthquakes occurring in close proximity to any occurred event in the record-breaking sequence.Incorporating this feature into the transition matrices results in modified transition matrices, as shown in Figs.2b  and 3b.We propose that accounting for the quiescence period in this manner opens additional options such as feedback effects on the non-linear dynamics of weakly coupled oscillators.

Synchronization
Scholz (2010) argued for the role of synchronization in fault interactions and earthquake clustering and for the usefulness of the Kuramoto model.Kuramoto (1975) proposed a mathematical model of phase oscillators interacting at arbitrary intrinsic frequencies and coupled through a sine of their phase differences.He suggested the following equations for each oscillator in the system:  where θ i is the phase of the ith oscillator, θi (t) is the first derivative of the phase of the ith oscillator with time, ω i is the natural frequency of the oscillator, K i is the strength of coupling of the ith oscillator to other oscillators and N is the size of the population of the oscillators.The frequencies ω i are chosen from a uniform distribution.Kuramoto (1975) demonstrated that synchronization was accomplished in the case of mean-field coupling with in the above equation.We can describe the Kuramoto model in a simpler form by introducing the complex-valued order parameter r(t): where (t) is the average phase and r(t) honours 0 ≤ r(t) ≤ 1.The expression of the Kuramoto model becomes The collective behavior of all the oscillators is monitored by examining the time evolution of the order parameter, r (Kuramoto, 1975;Strogatz, 2000;Pikovsky et al., 2003;Strogatz, 2003).The order parameter can assume values in the range 0 to 1 including the limits.From this, it is obvious that each oscillator is connected to the common average phase with the coupling strength is given by Kr.A value of "0" for r corresponds to total incoherence, i.e. no phase locking of the phases of the oscillators; a value of "1" for r corresponds to full coherence, i.e. phase locking of all the phases of the oscillators.The time evolution of the Kuramoto model can be monitored either by looking at the polar plots of the phases on a unit circle (Kuramoto, 1975) or by following the plot of the order parameter, r, as a function of the coupling strength, K. Acebrón et al. (2005) have provided a comprehensive review of the Kuramoto model.
For the stability of the solution from the Kuramoto model, use of a large population of oscillators for calculability in the thermodynamic limit is a pre-requisite.Over the last decade, efforts have gone into considering a finite number of oscillators satisfying the original conditions of the Kuramoto model.Easing the restrictions on the interaction model can be cast as an investigation of synchronization on complex networks.This would allow one to relate the complex topology and the heterogeneity of the network to the synchronization behavior.
We rewrite the original Kuramoto model for a complex network corresponding to undirected and directed graphs as where K ij is the coupling strength between pairs of connected oscillators and a ij refers to the elements of the adjacency or connectivity matrix.Much effort has gone into understanding the role of the coupling strength (Hong et al., 2002;Arenas et al., 2008;Dörfler et al., 2013) in the synchronization behavior of small-world and scale-free graphs.Here, we leave the coupling strength term a constant, unlike in the model under the thermodynamic limit where the size of the population, N, enters explicitly in the coupling strength term as a divisor.The structure of the adjacency matrix decides essentially the nature of the interaction term made up of the sine coupling of the phases.Vasudevan and Cavers (2014a) have investigated the synchronization behavior of the random graphs under different rewiring probabilities and the scalefree graphs from a spectral graph theory point of view.These studies did not include a study on the effect of clustering on the synchronization.In this regard, the work of McGraw and Menzinger (2005) is quite appealing.They conclude that for random networks and scale-free networks, increased clustering promotes the synchronization of the most connected nodes (hubs) even though it inhibits global synchronization.We see the role of the effect of clustering on the nature of synchronization behavior in earthquake sequencing studies and will conduct a separate study.Whether or not we reach similar conclusions for directed graphs, we have recently investigated synthetic networks that mimic real data structures (Vasudevan and Cavers, 2014a).In this regard, it is worth mentioning that synchronization of Kuramoto oscillators in directed networks has been subjected to a detailed study (Restrepo et al., 2006).

Chimera states
While the synchronization and asynchronization studies on earthquake sequencing are important in terms of the Kuramoto model given in Eq. ( 9), very little attention has been paid to the co-existence of synchronized and asynchronized states or the chimera states.Kuramoto and Battogtokh (2002) and Abrams and Strogatz (2004) paved the way for such a study by including the non-local effects in the Kuramoto model, as expressed in Eqs.(1) and (2).Non-local effects mean simply the inclusion of geometry effects.For the global seismicity map considered in this study, we generated the shortest-path distance matrices with and without the inclusion of recurrences (Figs. 2c,3c,2d,and 3d).The shortestpath algorithm encapsulates both the cascading effects of earthquakes and the negation of long-range distance effects.
In this study, we kept the global coupling strength constant and allowed the non-local coupling strength, κ, to vary from one simulation to the next one, similar to what was done in the recent work of Zhu et al. (2014).
Symmetry-breaking phenomena like chimera states have also been observed for two-cluster networks of oscillators with a Lorentizian frequency distribution (Montbrió et al., 2004) for all values of time delay.A crucial result by Laing (2009aLaing ( , 2009b) ) extends the previous observation to oscillators with heterogeneous frequencies.Also interesting to observe in this regard is that these heterogeneities can lead to new bifurcations, allowing for alternating synchrony between the distinct populations over time.Ko and Ermentrout (2008) demonstrated the presence of chimera-like states when the coupling strengths were heterogeneous.The last study used coupled Morris-Lecar oscillators.Although there is overwhelming evidence for the existence of chimera states in the presence of time delay or phase lag, all of our initial Kuramoto model simulations on the directed graph transition matrices and the associated shortest-path distance matrices included a constant phase lag only.
A postulation for the existence of evolving chimera states in data from earthquake catalogues has certain implications.For instance, it would pave the way to understanding the evolving alterations in stress-field fluctuations in fault zones frequented by earthquakes.Also, it would suggest a need to consider steps to quantify partially or fully the ratio of the number of synchronized oscillators to the total number of oscillators.The steps would involve extensive testing of the dependence of the parameters and additional mathematical models.We interpret the zones with synchronized oscillators as the ones being susceptible to earthquakes and the zones with asynchronized oscillators as the ones going through a quiescence period.The hope is that confirmation of chimera states in earthquake sequencing would signal a possible use for earthquake forecasting studies.

Simulation results and analysis
The Kuramoto model simulation with non-local coupling effects (κ = 0.10) with a phase lag, as expressed in Eq. ( 1) for a 128 × 128 grid transition probability, and the corresponding shortest-path distance matrices, lead to snapshots of three attributes: (i) the phase profile, (ii) the effective angular velocities of oscillators, and (iii) the fluctuation of the instantaneous angular velocity of oscillators.We did not sort the results according to an increase in the values of the attributes.Figure 4a-c show that, for a case of no recurrences, there exists a chimera state.The ensemble averages from the last 1000 time steps of the 200 000 time steps in the numerical simulations reveal the co-existence of synchronous and asynchronous oscillators.This means that some of the cells in the grid show a synchronous behavior and some others do not.In this particular case of no recurrences (Fig. 4), the number of synchronous oscillators to the number of asynchronous oscillators is large.In the case of recurrences, as shown in Fig. 5, this ratio is much larger.Also, the chimera pattern of the synchronized and asynchronized components of the oscillators is similar to what was observed by Abrams and Strogatz (2004).Figures 4a and 5a are the first evidence of the possible existence of a chimera state in earthquake sequencing.Figures 4b, c, 5b, and c confirm this.
Going from the 128 × 128 grid to the 1024 × 1024 grid, there are more non-zero cells with multiplicity of earthquakes at least 1.The number of oscillators is substantially larger, 1693 vs. 7697.Figures 6 and 7 reveal the chimera state as the steady-state solution to the non-linear dynamics on weakly coupled oscillators without and with recurrences for the 1024 × 1024 grid with the non-local coupling coefficient, κ = 0.1.Figure 8 shows the behavior of the chimera index as a function of the number of time steps for the nonrecurrent 128 × 128 grid with the non-local coupling coefficient, κ, set at 0.10.The purpose of this figure is to demonstrate that the asymptotic behavior of the chimera index with an increase in the number of time steps could be used to look at the steady-state solution of the Kuramoto model.
We looked at the influence of the non-local coupling coefficient, κ, on the ratio of the number of coherent oscillators to the total number of oscillators for both the 128 × 128 and 1024 × 1024 grids without recurrences in Fig. 9.A similar observation is made for the case of recurrences.As the nonlocal coupling coefficient, κ, increases from 0.01 to 1.0, the ratio decreases.For values of κ approaching 0, the non-local Kuramoto model acts as a simple Kuramoto model in that there is full synchronization for the global coupling parameter, A (or used as K in the literature), of 1.0.What is surprising to begin with is that, as κ approaches 1, the steady-state solution becomes more asynchronized.Investigations on the effect of the non-local coupling effect parameter, κ, on the steady-state solution of the phase angle distribution in the chimera state (Figs. 10a,b,11a and b) suggest that for both the 128 × 128 grid and the 1024 × 1024 grid, for larger κ values, the number of asynchronous oscillators is larger, and for smaller κ values, the presence of synchronous oscillators becomes dominant.For in-between values, i.e. between 1.0 and 0.03, the nature of the chimera states changes.
The outcome of each one of the simulations described for both the non-recurrence and recurrence cases contains synchronous and asynchronous vectors.Mapping these vectors on the respective grids (128 × 128 or 1024 × 1024 grids) should reveal the "non-readiness or readiness" cells or zones for earthquakes.One such map for a 128 × 128 grid without recurrences for κ = 0.10 is shown in Fig. 12.This qualitative description of the evolutionary dynamics of the earthquake sequencing is highly instructive.

Conclusions and future work
Earthquake sequencing is an intriguing research topic.The dynamics involved in the evolution of earthquake sequencing are complex.Very much has been understood, and yet the evolving picture is incomplete.In this regard, the work of Scholz (2010) acted as a catalyst in us investigating the synchronization aspect of earthquakes using the Kuramoto model.To name a few, the works of Vieira (1999), Rundle et al. (2002Rundle et al. ( , 2003)), Kuramoto and Battogtokh (2002), Abrams andStrogatz (2004), andLaing (2009a, b) have helped us take this step forward with this work.We summarize below the main points of this paper and also point out the direction in which we are going: 1. Earthquake sequencing from the IRIS earthquake catalogue browser can be expressed as a transition matrix of a directed graph.3.For a non-local coupling strength, κ, of 0.10, the Kuramoto model yields chimera states as a steady-state solution, i.e. co-existence of synchronized and asynchronized states.This is true for all the grid sizes considered.Differences exist in the ratio of the number of coherent oscillators to the number of incoherent oscillators.4. As the non-local coupling strength, κ, is lowered from 1.0 to 0.01, there is a general tendency towards an increase in synchronization, as is expected.While this general trend is observed for directed graphs generated from grids of orders 128, 192, 256, and 512, the graph from the 1024 × 1024 grid reveals the presence of the chimera state.
5. As the non-local coupling strength is increased from 0.1 to 1.0, there is a steady increase in the asynchronous behavior.
6.The recurrence results support the presence of chimera states for both 128 × 128 and 1024 × 1024 grids.However, it is quite intriguing to find out that the asynchronous oscillators come from a sub-set of the oscillators in both cases.
7. There is still a nagging question about which non-local coupling coefficient would be an ideal candidate for understanding the global stress-field fluctuations.Figure 12 illustrates an example of how a chimera state could be displayed on the map grid.Imposing geophysical and geodetic constraints on the earthquake zones in terms of heterogeneity of the natural frequencies would provide a quantitative answer to the above question.
1.In general, the hypothesis that all networks of earthquake faults around the globe go through full synchronization still needs to be strongly tested.On the other hand, the prevalence of chimera states or multi-chimera states is an attractive option to understand the earthquake sequencing.
2. We believe that there is, now, a mechanism available to us to explore and seek an answer to the non-linear dynamics of earthquake oscillations.
Needless to say, the role of the parameters such as the heterogeneity of the oscillators as expressed in the natural frequency of the oscillators, the variability of the time-delay corrections instead of a constant time delay, and the heterogeneity of the non-local coupling strength and the global coupling strength in the present Kuramoto model, remain to be investigated.Work is currently in progress.

Figure 1 .
Figure 1.Partitioning of the global seismicity map: (a) 128 × 128 gridding of the latitude-longitude map.(b) 1024 × 1024 gridding of the latitude-longitude map.Earthquakes of magnitudes exceeding or equal to M w = 5.5 and location depth not exceeding 70 km for the time period from January 1970 to September 2014 constitute the glacial seismicity map.Earthquake information was downloaded from IRIS (Incorporated Research Institutions for Seismology).The earthquake frequency used in the maps is plotted on a log(log) display scale, with larger circles representing higher frequencies.

Figure 2 .
Figure 2. 128 × 128 gridded map: (a) transition probability matrix without recurrences.(b) Transition probability matrix with recurrences.(c) Shortest-path distance matrix without recurrences.(d) Shortest-path distance matrix with recurrences.In (a) and (b), the transition frequencies used in the maps are plotted using a log(log) display scale, with larger circles representing higher frequencies.

Figure 3 .
Figure 3.A 1024 × 1024 gridded map: (a) transition probability matrix without recurrences.(b) Transition probability matrix with recurrences.(c) Shortest-path distance matrix without recurrences.(d) Shortest-path distance matrix with recurrences.In (a) and (b), the transition frequencies used in the maps are plotted using a log(log) display scale, with larger circles representing higher frequencies.

Figure 4 .
Figure 4. Three attributes of a chimera state of the 1693 oscillators for a 128 × 128 gridded map without recurrences using κ = 0.10.(a) Stationary phase angle.(b) Effective angular velocity.(c) Fluctuations in instantaneous angular velocity.

Figure 5 .
Figure 5. Three attributes of a chimera state of the 1693 oscillators for a 128 × 128 gridded map with recurrences using κ = 0.10.(a) Stationary phase angle.(b) Effective angular velocity.(c) Fluctuations in instantaneous angular velocity.

Figure 6 .
Figure 6.Three attributes of a chimera state of the 7697 oscillators for a 1024 × 1024 gridded map without recurrences using κ = 0.10.(a) Stationary phase angle.(b) Effective angular velocity.(c) Fluctuations in instantaneous angular velocity.

Figure 7 .
Figure 7. Three attributes of a chimera state of the 7697 oscillators for a 1024 × 1024 gridded map with recurrences using κ = 0.10.(a) Stationary phase angle.(b) Effective angular velocity.(c) Fluctuations in instantaneous angular velocity.

Figure 8 .
Figure 8. Chimera index as a function of time steps for the 128 × 128 grid without recurrences for κ = 0.10.

Figure 9 .
Figure 9. Influence of the non-local coupling coefficient parameter, κ, on the ratio of the number of synchronized oscillators to the total number of oscillators for both the 128 × 128 and the 1024 × 1024 grids without recurrences.

Figure 12 .
Figure 12.Chimera-state map of the synchronous and asynchronous oscillators as a steady-state solution for a non-recurrence case.The non-local coupling coefficient parameter, κ, is 0.1.Blue dots refer to the asynchronous oscillators and red dots to the synchronous oscillators.

Table 1 .
Grid sizes and the number of oscillators corresponding to non-zero cells.