Microscopic Cross-Correlations in the Finite-Size Kuramoto Model of Coupled Oscillators
aa r X i v : . [ n li n . AO ] J u l Microscopic Cross-Correlations in the Finite-Size Kuramoto Model of CoupledOscillators
Franziska Peter, Chen Chris Gong, and Arkady Pikovsky Institute of Physics and Astronomy, University of Potsdam,Karl-Liebknecht-Straße 24-25, 14476 Potsdam, Germany Department of Control Theory, Nizhny Novgorod State University,Gagarin Av. 23, 606950, Nizhny Novgorod, Russia (Dated: July 29, 2019)Super-critical Kuramoto oscillators with distributed frequencies separate into two disjoint groups:an ordered one locked to the mean field, and a disordered one consisting of effectively decoupledoscillators – at least so in the thermodynamic limit. In finite ensembles, in contrast, such clearseparation fails: The mean field fluctuates due to finite-size effects and thereby induces order inthe disordered group. To our best knowledge, this publication is the first to reveal such an effect,similar to noise-induced synchronization, in a purely deterministic system. We start by modelingthe situation as a stationary mean field with additional white noise acting on a pair of unlockedKuramoto oscillators. An analytical expression shows that the cross-correlation between the twoincreases with decreasing ratio of natural frequency difference and noise intensity. In a deterministicfinite Kuramoto model, the strength of the mean field fluctuations is inextricably linked to thetypical natural frequency difference. Therefore, we let a fluctuating mean field, generated by afinite ensemble of active oscillators, act on pairs of passive oscillators with a microscopic naturalfrequency difference between which we then measure the cross-correlation, at both super- and sub-critical coupling.
I. INTRODUCTION
Synchronization – the mutual adjustment of frequen-cies among weakly coupled self-sustained oscillators – isa prominent example of the emergence of order in out-of-equilibrium systems in physics, engineering, biology,and other fields [1–3]. In large ensembles, it appears asa non-equilibrium phase transition, where the organizingaction of sufficiently strong mutual coupling wins overthe disorganizing action of to the diversity in natural fre-quencies. The paradigmatic model of this phenomenon,created by Kuramoto [1, 4], is fully solvable in the ther-modynamic limit [1, 5]. The characteristic feature of theKuramoto-type synchronization transition is the coex-istence of two subgroups of oscillators in the partiallysynchronized state: the oscillators in the ordered groupare locked by the mean field and coherently contributeto it, while the disordered units are not locked and ro-tate incoherently. With increasing coupling strength, theformer group grows in size, as more and more oscillatorsare locked by the increasing mean field.The qualitative features, established in the thermody-namic limit, remain approximately valid for finite ensem-bles. Here, similar to finite-size effects in equilibriumphase transitions, the order parameter, i. e. the macro-scopic mean field, fluctuates with an amplitude that de-pends on the ensemble size in a nontrivial way [6–10].These fluctuations are most pronounced close to the crit-icality, and can be attributed to weak chaoticity of thefinite population dynamics [11, 12].The goal of this paper is to show that the finite-sizefluctuations of the mean field have an additional effecton the population – quite counter-intuitively an order-ing effect: the disordered oscillators become correlated pairwise, while in the thermodynamic limit the cross-correlations disappear. Below, we consider two basic se-tups to show this phenomenon. First, we study a popula-tion in the thermodynamic limit, but with the mean fieldbeing subject to external white noise fluctuations (Sec-tion II). This ideal setup allows for an analytic solution,showing the dependence of the cross-correlation betweenthe oscillators on the fluctuation intensity and on thenatural frequency difference. In the second setup (Sec-tion III), we numerically quantify the cross-correlationdue to the intrinsic finite-size-induced fluctuations of themean field, first for super- than for sub-critical coupling.This latter case is similar to other organizing macroscopicmanifestations of finite-size fluctuations such as finite-size-induced phase transitions [13, 14] and stochastic res-onance [15]. In fact, this ordering action of finite-sizefluctuations can be qualitatively traced to the effect ofsynchronization by common noise, known for identicaland nonidentical oscillators, which are either coupled oruncoupled [16–19].
II. MEAN FIELD WITH EXTERNALFLUCTUATIONS IN THERMODYNAMIC LIMITA. Stationary mean field
Before discussing the mean-field model with externalfluctuations, we first quantify ordered and disorderedstates in the Kuramoto model of mean-field coupled os-cillators in the thermodynamic limit where no fluctu-ations are present. The model is formulated as fol-lows: Oscillators are described by their phases ϕ , andare coupled via the complex mean field Z ≡ R e i Φ = R π dϕ R ∞−∞ d Ω P ( ϕ | Ω) g (Ω)e iϕ as˙ ϕ = Ω + ε Im( Z e − iϕ ) . (1)Here Ω are natural frequencies distributed according toa unimodal density g (Ω), and P ( ϕ | Ω) is the probabilitydensity of oscillators with natural frequency Ω.The theory of synchronization, developed by Ku-ramoto [1], predicts the existence of a critical value ofthe coupling constant ε , beyond which the macroscopicmean field and the frequency of the global phase assumeconstant values, R > | Ω − Ω | < εR are locked by the mean field, i. e. rotate with ¯Ω, andconstitute the ordered part of the population. Oscil-lators at the tails of the natural frequency distributionwith | Ω − Ω | > εR each rotate with a different averagefrequency and constitute the disordered part.To quantify order and disorder in the system, we cal-culate the pairwise cross-correlations between the oscilla-tors. First, we perform a shift to the mean field referenceframe ˜ ϕ := ϕ − Φ. In the ordered part, the oscillatorshave constant phases ˜ ϕ and are thus perfectly correlated.To calculate cross-correlations in the disordered part, thephases ˜ ϕ cannot be used directly, because their probabil-ity density on the circle is not uniform: it is proportionalto ˜ P ( ˜ ϕ | Ω) ∼ | Ω − Ω − εR sin ˜ ϕ | − [3], i.e., is a wrappedCauchy distribution. This distribution is fully character-ized by the first harmonic z = h e i ˜ ϕ i = iq (1 − p − q − ),where q = (Ω − Ω) / ( εR ) and h·i is the average over os-cillator phases ˜ ϕ with density ˜ P ( ˜ ϕ | Ω). The phases ˜ ϕ canbe transformed to uniformly distributed phase variables ψ by virtue of a M¨obius transform [20]exp[ iψ ] = (exp[ i ˜ ϕ ] − z ) (cid:14) (1 − z ∗ exp[ i ˜ ϕ ]) . (2)Straightforward calculations show that ˙ ψ = ν = (Ω − Ω) (cid:0) − q − (cid:1) / , where ν denotes the observed frequencyof the oscillator. Because the transformed phases ψ ro-tate uniformly, with their frequency ν now only depend-ing on intrinsic frequency Ω, we can straightforwardly ap-ply the synchronization index – a measure for the cross-correlation of two phases [21] – as γ = |h exp[ i ( ψ − ψ )] i| . (3)For two phase variables uniformly rotating with differentfrequencies ν and ν , this yields zero, thus the pairwisecross-correlation vanishes exactly in the disordered do-main, independently of how close the parameters of theoscillators are. Conversely, this cross-correlation is ex-actly 1 in the order domain. B. Mean field with external fluctuations
Our next goal is to show that non-vanishing cross-correlations appear in the disordered part if the mean field contains added external noise. We consider a sim-ple modification of the Kuramoto model (1): Z = R exp[ i Ω t ] + ˜ σ [ ξ ( t ) + iξ ( t )], with constants R, Ω andGaussian random processes ξ ( t ) , ξ ( t ) of noise strength˜ σ with h ξ i ( t ) ξ j ( t ′ ) i = 2 δ ij δ ( t − t ′ ). We assume noise to beweak, so that it does not significantly change the individ-ual statistical properties of the oscillators (the distribu-tion remains a wrapped Cauchy distribution, see Ref. [22]for a quantification of small deviations due to weak Gaus-sian noise). However, as we will see, it induces cross-correlations in the disordered region. Performing thesame transformation to obtain uniformly rotating phasevariables ψ as above, we obtain for the transformed phasevariables ψ a set of Langevin equations with commonnoise terms η , η :˙ ψ = ν − ( a + b sin ψ ) η ( t ) + c cos ψ η ( t ) . (4)Here η ( t ) , η ( t ) are mutually uncorrelated Gaussianwhite noise forces common to all transformed phases ψ , h η i ( t ) η j ( t ′ ) i = 2 δ ij δ ( t − t ′ ); ν = (Ω − Ω) (cid:0) − q − (cid:1) / is the observed frequency as above; parameters a = ε R ˜ σ/ν , b = ε ˜ σ (Ω − Ω) /ν , and c = ε ˜ σ are the effectivenoise strengths.We now consider two oscillators ψ , ψ from the set (4).We assume parameters of these oscillators to be close: ν = ν + ρ/ ν = ν − ρ/
2, with ρ ≪ ν ; and similarfor a , and b , . To calculate the cross-correlation func-tion (3), we first write, starting from (4), the Langevinequations for the difference α = ψ − ψ and the sum β = ψ + ψ of the transformed phase variables:˙ α = ρ − ( a − a ) η ( t ) − ( b + b ) sin α β η ( t ) − ( b − b ) cos α β η ( t ) − c sin α β η ( t ) , ˙ β =2 ν − ( a + a ) η ( t ) − ( b + b ) cos α β η ( t ) − ( b − b ) sin α β η ( t ) + 2 c cos α β η ( t ) . This system yields a Fokker-Plank equation for thedensity W ( α, β, t ), which, by virtue of averaging overthe fast rotating variable β with the method of multi-ple scales [23], can be reduced to the following equationfor the density of the phase difference w ( α, t ): ∂∂t w + ρ ∂∂α w − (cid:18) ( b + b ) c (cid:19) ∂ ∂α [(1 − cos α ) w ] =( b − b ) ∂ ∂α [(1 + cos α ) w ] + ( a − a ) ∂ ∂α w . (5)The terms on the r.h.s. of this equation are of secondorder in the small parameter ρ , and therefore we canneglect them. As a result, only the weighted sum of noiseterms σ := ( b + b ) / c ≈ ε ˜ σ { − Ω) /ν ] } is relevant.The stationary solution of this equation can bestraightforwardly written as an integral; the cal-culation of the cross-correlation γ = (cid:12)(cid:12) h e iα i (cid:12)(cid:12) = | R π − π w ( α ) exp( iα ) d α | reduces to a nontrivial integra-tion, which nevertheless can be expressed explicitly: γ = 1 + 4 d [ci (2 d ) + si (2 d )] − d [ci(2 d ) sin(2 d ) − si(2 d ) cos(2 d )] , (6)where ci and si are the cosine and sine integral func-tions, and the ratio d between the frequency mismatch ρ and the noise strength σ , d = ρ/σ , is the only param-eter. This cross-correlation function (cf. Fig. 4, blacksolid line) tends to 1 for d ≪ γ ∼ /d as d ≫
1. Thus, our main analytical result (6) showsthat common external noise added to the mean field in-duces cross-correlations in the disordered domain, witha characteristic cross-correlation length proportional tothe noise intensity ρ ∼ σ .The physical explanation of this cross-correlation liesin the stabilizing effect of the common noise: its actionon an oscillator leads to a negative Lyapunov exponent,which results in complete synchronization of identical os-cillators [16, 17]. For nonidentical oscillators, the differ-ence in the natural frequencies prevents complete syn-chrony, but the phases are most of the time kept close toeach other by noise, with occasional fast phase slips [24]that account for the observed frequency difference. III. FINITE-SIZE MEAN FIELD FLUCTUATION
As demonstrated above, micro-scale cross-correlationsappear in the disordered domain of the Kuramoto modelin the thermodynamic limit with external mean fieldnoise. A natural question arises, if also the intrinsic orderparameter fluctuations in deterministic finite ensemblesgenerate such micro-scale cross-correlations. The essen-tial parameter here is ensemble size N . Simple estima-tions based on the theory above show that the cross-correlations are rather small between typical pairs os-cillators: For an ensemble of size N , the characteristicfrequency mismatch between the oscillators is ρ ∼ N − .However, in order to create sizeable fluctuations of themean field, N must be small. If one assumes σ ∼ N − ,then a typical value of the parameter d = ρ/σ will beclose to unity, which is too large for the cross-correlationsto be observable (see Fig. 4). A. A model with active and passive oscillators
The size of the ensemble dictates not only the size ofthe intrinsic fluctuations of the mean field, which tends tohave an organizing effect on the phases, it also determinesthe typical natural frequency mismatch of a given pair ofoscillators, which tends to have an disorganizing effecton their phases. According to (6), the cross-correlation between the phases with typical natural frequency dif-ference for typical mean field fluctuation intensity (bothdepending on N ) are small. The probability to find apair with much smaller than average natural frequencydifference is high, but then again it is difficult to disen-tangle different effects on such singular pairs. However,this problem can be resolved if the fluctuation level (orthe effective noise strength) is decoupled from the rangeof natural frequency differences.To decouple the two opposing effects, we introduce amodification of the Kuramoto model, where the oscilla-tors are of twoc) types: active ones φ j ( j = 1 , , . . . , N )with natural frequencies Ω Aj , and passive ones (tracers) ϕ k ( k = 1 , , . . . , M ) with natural frequencies Ω k . Theoscillators of both types obey the same equation (1).However, only the active oscillators contribute to themean field: Z = R exp[ i Φ] = N − P Nj =1 exp[ iφ j ]. Here N is the number of active oscillators, while the numberof passive ones M can be arbitrary (and they can haveany distribution of frequencies). One can say that passiveoscillators “test” the mean field created by active oscil-lators, similar to how ideal fluid tracers “test” the flowof a fluid. The passive oscillators do so at different fre-quencies, especially at those not presented in the activeset. Similar technique has been used in [25] to determinethe frequency of chaotic signals via locking.Equivalently, the system of active-passive oscillatorscan be considered as a large network˙ ϕ k = Ω k + εN X j K kj sin( ϕ j − ϕ k ) , (7)where K kj = 1 if the phase ϕ j belongs to the activeset, and K kj = 0 otherwise. Experimentally, such a cou-pling has been directly implemented in a set of 2816 op-tically coupled periodic chemical Belousov-Zhabotinskyreactors [26].Similar setups are often used in systems with long-range interactions, for example, in restricted N -bodyproblems in gravitational systems. Heavy bodies suchas planets, stars and galaxies contribute to the gravita-tional field in which they move, while other, lighter parti-cles move in the same field, but their contribution to thefield is negligible. In a more general context of interdisci-plinary applications of complex systems, the division intoactive and passive agents occurs by itself in macro-socialopinion formation processes. In social media, a few for-ward thinkers (or influencers) lead the public discourseby writing texts and comments, while the opinions of alarge number of passive users (followers) remain hidden:they follow the discussion without contributing to it [27]. B. Fluctuations beyond the synchronizationtransition
As we show below, using tracers, the micro-scale cross-correlations (and other interesting features) and can be ϕ k a) ϕ k b) Ω k −3−2−10123 ϕ k oldnew c) FIG. 1. Snapshot of passive oscillators (black dots) and ac-tive oscillators (large red dots) phases plotted against naturalfrequencies for a finite ensemble of 50 active oscillators (notall shown). The mean field is created by the active oscillatorswith natural frequencies drawn randomly from a Gaussian dis-tribution (same sample as in Fig. 2) at slightly super-criticalcoupling ε = 1 .
85 (same in Fig. 2). Three zoom levels of fac-tor 10 are marked by shaded areas. Arrows in c) point atΩ k ≈ . k ≈ .
72 for a new and an old phase slips,respectively. See Supplementary Material at [URL will beinserted by publisher] for an animated version of this figure.One second in the video equals one time unit. easily detected. In this subsection, we illustrate thesefeatures for a partially synchronized state of the finite-size Kuramoto model, and in subsection III C for a statebelow the synchronization transition.First, we give a qualitative picture of the cross-correlations. Fig. 1 shows a snapshot of phases for apopulation of N = 50 active oscillators together witha set of M = 5 · tracers. The natural frequenciesof the active ones are sampled from a normal distribu-tion with zero mean and unit variance, for which thecritical coupling constant in the thermodynamic limit is ε c ≈ .
6. The micro-scale ordering effect becomes evi-dent if one zooms in to increasingly smaller scales, frompanel a) to panel c). Panel c) shows the characteristiccorrelated state profile of the tracers’ phases, consistingof ruptured nearly horizontal bars. A bar is formed dueto the ordering action of the fluctuations of the meanfield, which synchronize passive oscillators with close fre- quencies. Ruptures appear when oscillators with higherfrequencies make an additional rotation (a phase slip)with respect to oscillators with smaller (but similar) fre-quencies. In Fig. 1c) one can clearly see a fresh phase sliparound Ω ≈ . ≈ . k & .
0) (Fig. 1a). Here thetracers are synchronized to the active units (like the satel-lites are trapped by their planet’s gravitational field).This is because the fluctuations of the mean field are infact not completely random, but contain relatively strongnearly periodic components from the non-entrained ac-tive oscillators. These components suffice, at least forsmall ensembles, to fully entrain tracers with natural fre-quencies close to a common active oscillator.Now we quantify the cross-correlations illustrated inFig. 1. We calculate the cross-correlation coefficient γ (Ω , ∆) for two tracers with natural frequencies Ω − ∆ / , Ω + ∆ / ψ . First, we calculatethe time dependent difference between the tracer phasesand the mean field phase ˜ ϕ ( t ) = ϕ ( t ) − Φ( t ). We then av-erage these phases over time, z = h exp[ i ˜ ϕ ] i , which givesthe empirical value of the parameter characterizing thewrapped Cauchy distribution of ˜ ϕ . Then, the M¨obiustransform (2) is applied. To check that we indeed ob-tained the uniformly distributed phase variable ψ , wecalculate the first harmonics z ψ = h exp[ iψ ] i and compareit to z (Fig. 2a): one can see that indeed the transforma-tion yields a uniformly distributed set of phase variables(within reasonable tolerance), because the amplitudes ofthe time averages of their first harmonics become veryclose to zero after the transformation.In Fig. 2 b) we show the observed frequencies of thetracers and the active oscillators. One can clearly seesynchronized neighborhoods of active units as plateausin this graph. Outside the plateaus, the tracers arenot locked and their observed frequency varies contin-uously with their natural one. It is in these domains out-side the plateaus, that the micro-scale cross-correlationscan be observed and measured, as illustrated in Fig. 2c). Here we show values of the cross-correlation coef-ficient γ (Ω , ∆) for several values of frequency mismatch∆: cross-correlation is nearly perfect for ∆ . − , whilefor ∆ & .
01 the values of the coefficient typically do notexceed 0 . |z k | z ˜ϕk z ψk o b s e r v e d f r e q u e n c i e s passive (cid:0) ˙ϕ k (cid:3) active (cid:1) ˙φ j (cid:0) natural frequencies Ω k , γ Ω Aj ∆ −6 −5 −4 −3 −2 a)b)c)FIG. 2. Observed cross-correlation as defined by Eq. (3) ofpassive oscillators coupled to the mean field Z of 50 active os-cillators against their natural frequencies (mean h R i = 0 . σ R = 3 . · − ), shown for various levels of frequencymismatch ∆ between passive oscillators. a): Absolute valueof the (time-averaged) first harmonic of individual passive os-cillators before ( z , blue circles) and after ( z ψ , red squares)the M¨obius transform. z ψ drops to near zero after the trans-formation, which shows that ψ is now uniformly rotating. b):Observed vs. natural frequencies of passive (line) and activeoscillators (dots). Gray vertical lines in panels b), c) mark thenatural frequencies of the active oscillators. c): Synchroniza-tion index γ between pairs of passive oscillators with naturalfrequencies Ω k ± ∆ / C. Fluctuations below the synchronizationtransition
To show that the effect of microscopic cross-correlations occurs for subcritical values of coupling con-stant ε as well, we present the calculations of the cross-correlations for an ensemble of N = 50 oscillators withnatural frequencies randomly sampled from the standardnormal distribution, and ε = 1 in Fig. 3. For this rela-tively small coupling, the complex mean field fluctuatesaround zero, see e. g. [10]. Therefore, a transformationof the phases to uniformly rotating ones ( ˜ ϕ → ψ ) is un-necessary, contrary to the case of a non-vanishing meanfield at stronger coupling (Figs. 1, 2). For better visi- γ −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 natural frequencies Ω k , -6 -5 -4 -3 -2 -1 l og ( − γ ) Ω Aj a)b)FIG. 3. The cross-correlations of passive oscillators coupled toa mean field of 50 active ones at ε = 1 (same natural frequencysample as Fig. 2). The mean field Z fluctuates around zero,see [10]. Panel a): cross-correlations in the linear scale; panelb): the same in logarithmic scale to resolve the region γ . bility of both high and low cross-correlations, we presentin Fig. 3 the cross-correlation constant γ (Ω , ∆) in lin-ear and logarithmic scales. The multiple locked regionsrelate to the frequencies of active oscillators. The cen-tral region around Ω ≈ S between the peak value and the next largestmaximum of the cross-correlations. This quantity de-creases with the system size of the active ensemble. Fur-thermore, we have chosen only non-locked passive units.One can see that the scaling relation γ = γ (∆ /S ) followsat least qualitatively the theoretical curve (6), althougha huge diversity of the observed cross-correlations, dueto the “coloredness” of the intrinsic finite-size mean fieldfluctuations, is also evident. -4 -3 -2 -1 d -1 γ N osc ε =12550100200 ε =1.85 2550100200 PSfrag replacementslog ( d ) γ ( d )FIG. 4. Comparison of the observed cross-correlation γ insub- and super-critical ensembles ( ε = 1 and 1 .
85 with openand filled symbols, respectively) of different sizes N osc =25 , , ,
200 with the theoretical expression (6) (solid line).The data points were generated in experiments as in Figs. 2, 3.The median of all γ ’s for which γ (10 − ) < . d = ∆ /σ is determined from the average over σ = 2 ε D { | Ω k − ¯Ω | / ∆) } with diffusion coefficient D that was integrated from the auto-correlation function of R .For each ensemble size, only one set of natural frequencies isrepresented. IV. CONCLUSION
In summary, we have shown that the fluctuations ofthe mean field in the Kuramoto model, either externallyimposed on the ensemble of infinite size, or naturally in-duced in the finite-sized model, lead to the appearance ofcross-correlations in the disordered part of oscillator pop-ulations. These cross-correlations result from the com-petition between synchronization by common noise anddesynchronization due to the parameter differences (usu-ally the differences in the natural frequencies). We havedeveloped an analytical theory of these cross-correlationsfor a mean field being a constant (in a properly rotatingframe) plus additional Gaussian white noise, summarizedin expression (6). This theoretical result is directly ap-plicable to models similar to the Kuramoto model, e.g.,to the Kuramoto-Sakaguchi model, where the mean fieldof a population is subject to external fluctuations. In thederivation of (6) we explicitly restrict the mean-field cou- pling to the first harmonics of the oscillator phases only.The case of a more general Daido-type coupling functionrequires extra analysis, although qualitative argumentsimply that the cross-correlations will be observed thereas well.Furthermore, we have numerically characterized pair-wise cross-correlations between passive oscillators drivenby the intrinsically fluctuating mean fields in a finite en-semble at two different coupling strengths, super- andsub-critical, respectively. In both cases, the mean fieldcontains nearly periodic components, therefore there ex-ist locked regions where high values of pairwise cross-correlations arise due to resonant locking. Between theselocked regions, the cross-correlations decay with increas-ing frequency mismatch, which is in a qualitative agree-ment with the theory based on white noise. The rathergood quantitative agreement between white noise theoryand finite Kuramoto model surprises, as, due to the dom-inant nearly periodic components in the fluctuating meanfield, fluctuations in the finite Kuramoto model are farfrom being “white”.Overall, the effect is expected to be most pronounced insituations where finite-size fluctuations are anomalouslylarge (populations with equidistant natural frequenciesat sub-critical coupling appear, as our preliminary calcu-lations show, to belong to this class).We expect that this phenomenon is not restricted tothe mean-field coupling, and can be observed in otherlarge systems where synchronized and disordered sub-populations coexist. A prominent example here is achimera state in a one- or two-dimensional oscillatorymedium with long-range interactions [28]. A popula-tion of oscillators driven by two mean fields [29] alsodemonstrated nontrivial regimes with a coexistence ofordered and disordered subpopulations. Indeed, chimerastates are defined as the coexistence of coherent and non-coherent domains among identical oscillators, and finite-size fluctuations [30] may lead to cross-correlations in thedisordered domain; this issue is currently under consid-eration.
ACKNOWLEDGMENTS
This paper was developed within the scope of theIRTG 1740 / TRP 2015/50122-0, funded by the DFG/FAPESP. A. P. was supported by the Russian ScienceFoundation (grant No. 17-12-01534). The authors thankD. Goldobin, R. Cestnik, M. Wolfrum, O. Omel’chenko,D. Paz´o, H. Engel, and R. Toenjes for helpful discussions. [1] Y. Kuramoto,
Chemical Oscillations, Waves, and Turbu-lence , Springer Series in Synergetics ed., Vol. 19 (SpringerBerlin Heidelberg, 1984).[2] S. H. Strogatz,
Sync: The Emerging Science of Sponta-neous Order (Hyperion, NY, 2003). [3] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchroniza-tion: A Universal Concept in Nonlinear Sciences (Cam-bridge University Press, 2001).[4] J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort,and R. Spigler, Rev. Mod. Phys. , 137 (2005). [5] E. Ott and T. M. Antonsen, Chaos (2008).[6] H. Daido, Journal of Physics A: Mathematical and Gen-eral , L629 (1987).[7] H. Hong, H. Chat´e, H. Park, and L. H. Tang, Phys. Rev.Lett. , 1 (2007).[8] I. Nishikawa, K. Iwayama, G. Tanaka, T. Horita, andK. Aihara, Progress of Theoretical and ExperimentalPhysics , 1 (2014).[9] H. Hong, H. Chat´e, L.-H. Tang, and H. Park, Phys. Rev.E , 022122 (2015).[10] F. Peter and A. Pikovsky, Phys. Rev. E , 032310(2018).[11] A. Pikovsky and A. Politi, Lyapunov Exponents: A Toolto Explore Complex Dynamics (Cambridge UniversityPress, 2016) Chap. 10.[12] B. Gilad,
Synchronization of Network Coupled Chaoticand Oscillatory Dynamical Systems , Ph.D. thesis, Uni-versity of Maryland (2013).[13] A. S. Pikovsky, K. Rateitschak, and J. Kurths, Z. PhysikB , 541 (1994).[14] M. Komarov and A. Pikovsky, Phys. Rev. E , 020901(2015).[15] A. Pikovsky, A. Zaikin, and M. A. de la Casa, Phys.Rev. Lett. , 050601 (2002).[16] D. S. Goldobin and A. S. Pikovsky, Radiophysics andQuantum Electronics , 910 (2004).[17] D. S. Goldobin and A. Pikovsky, Phys. Rev. E ,045201(R) (2005).[18] K. H. Nagai and H. Kori, Phys. Rev. E , 065202 (2010). [19] A. V. Pimenova, D. S. Goldobin, M. Rosenblum, andA. Pikovsky, Scientific Reports , 38518 (2016).[20] The transformation ˜ ϕ → ψ is an example of the pro-tophase to phase transformation used in the data analysisof oscillatory systems [31].[21] K. Mardia and P. Jupp, Directional Statistics , Wiley Se-ries in Probability and Statistics (Wiley, 2009).[22] I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, andA. Pikovsky, Phys. Rev. Lett. , 264101 (2018).[23] A. Nayfeh,
Introduction to Perturbation Techniques , Wi-ley classics library (Wiley, 1981).[24] D. S. Goldobin, A. V. Pimenova, M. Rosenblum, andA. Pikovsky, The Eur. Phys. J. Spec. Topics , 1921(2017).[25] M. Rosenblum, A. Pikovsky, J. Kurths, G. Osipov,I. Kiss, and J. Hudson, Phys. Rev. Lett. , 264102(2002).[26] J. F. Totz, J. Rode, M. R. Tinsley, K. Showalter, andH. Engel, Nature Physics , 282 (2017).[27] J. Gerson, A. C. Plagnol, and P. J. Corr, Personalityand Individual Differences , 81 (2017).[28] Y. Kuramoto and D. Battogtokh, Nonlinear Phenom.Complex Syst. , 380 (2002).[29] X. Zhang, A. Pikovsky, and Z. Liu, Scientific Reports ,2104 (2017).[30] M. Wolfrum and O. E. Omel’chenko, Phys. Rev. E ,015201 (2011).[31] B. Kralemann, L. Cimponeriu, M. Rosenblum,A. Pikovsky, and R. Mrowka, Phys. Rev. E77