High-Order Phase Reduction for Coupled Oscillators
Erik Genge, Erik Teichmann, Michael Rosenblum, Arkady Pikovsky
HHigh-Order Phase Reduction for CoupledOscillators
Erik Gengel , Erik Teichmann , Michael Rosenblum , ArkadyPikovsky , Institute of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str.24/25, 14476 Potsdam-Golm, Germany Control Theory Department, Institute of Information Technologies, Mathematicsand Mechanics, Lobachevsky University, Nizhny Novgorod, Russia
Abstract.
We explore the phase reduction in networks of coupled oscillators in thehigher orders of the coupling parameter. For coupled Stuart-Landau oscillators, wherethe phase can be introduced explicitly, we develop an analytic perturbation procedureto allow for the obtaining of the higher-order approximation explicitly. We demonstratethis by deriving the second-order phase equations for a network of three Stuart-Landauoscillators. For systems where explicit expressions of the phase are not available, wepresent a numerical procedure that constructs the phase dynamics equations for asmall network of coupled units. We apply this approach to a network of three vander Pol oscillators and reveal components in the coupling with different scaling in theinteraction strength. a r X i v : . [ n li n . AO ] J u l igh-Order Phase Reduction for Coupled Oscillators
1. Introduction
Networks of coupled self-sustained oscillators are widely used to describe complexrhythmical systems in physics [1, 2], biology [3, 4], and other fields [5, 6]. In particular,such models are relevant for the description of laser [7] and nanomechanical [8] oscillatorarrays, coupled Josephson junctions [9] and spin-torque oscillators [10], power grids [11],the activity of neuronal populations [12], the interaction of different organs within ahuman body [13], cell assemblies [14], etc.One of the most famous theoretical tools for the analysis of coupled oscillatorsis the phase reduction [15–18]. This approach provides a recipe for a description ofcomplex high-dimensional oscillators by a single cyclic variable, phase ϕ , so that thedynamics of a network of N generally high-dimensional elements reduces to a set of N coupled one-dimensional differential equations for the phases. This reduction is basedon the slaving principle (the corresponding notion in mathematical literature is thenormal hyperbolicity): While the phases correspond to the neutral directions with zeroLyapunov exponents, the “amplitudes” (i.e. all other variables except for the phases)are stable and thus follow the phase dynamics. As a result, the full system’s dynamicsreduces to that on the N -dimensional torus spanned by the phases.This reduction allows studying general behaviours of coupled oscillators throughphase dynamics models, a prominent example here is the analytically tractableKuramoto model of all-to-all interconnected units, derived for a population of weaklycoupled oscillators close to the Hopf bifurcation point [15, 19, 20]. Generally, phasedynamics models are expected to be valid for arbitrary oscillators and weak tomoderate coupling. However, the existing theory employs a perturbative approach andprovides phase equations only in the first-order approximation in the coupling strength.Derivation of high-order corrections would extend the validity of the approach beyondthe weak-coupling limit and, thus, essentially increase our ability to analyse complexnetworks. Yet, for the moment it remains a theoretical challenge, in spite of numerousattempts [21–26].In addition to obvious advantages for theoretical studies, phase reduction providesa framework for model reconstruction from measurements and, hence, plays an essentialrole in experimental research. Reconstruction of the phase dynamics equations fromdata is much more simple than recovery of equations for the state variables because (i)phase equations are low-dimensional and (ii) have a simple universal form because theirright-hand sides are 2 π -periodic functions of the phases. Though precise estimationof the phases from scalar signals remains a topic of ongoing studies [27–29], phasemodels have been successfully recovered for several laboratory experiments as well asfor physiological systems [30–38]. In particular, this approach yields a way to tackle theconnectivity problem, i.e. to recover directional causal links between oscillatory sourcessolely from multivariate data.However, for more than two oscillators and moderate coupling the connectivity, as itappears in phase equations, generally differs from the connectivity as defined by physical igh-Order Phase Reduction for Coupled Oscillators
2. Theoretical analysis
We first briefly recall the basic ideas of the phase reduction in the first order in thecoupling parameter ε (see Refs. [16, 18] for more details). One starts by considering anoscillatory dynamical system d y dt = f ( y )possessing a stable limit cycle Y ( t ) = Y ( t + T ) in its state space. As is well-known, onthe limit cycle and in its basin of attraction one can define the phase ϕ = Φ( y ) whichgrows uniformly in time˙ ϕ = ω = ∂ Φ ∂ y f ( y )with basic frequency ω = 2 π/T . Notice, that on the limit cycle the system’s state isuniquely defined by the phase: y = Y ( ϕ ). Outside of the limit cycle, this is not true:one also has to know the deviation from the cycle.In the context of oscillatory networks, one considers many coupled oscillators thatwe label by index k : d y k dt = f k ( y k ) + ε G k ( y , y , . . . ) . (1) igh-Order Phase Reduction for Coupled Oscillators ε is the small parameter governing the strength of the coupling. The equation forthe phases is obtained by exploiting their definition ϕ k = Φ k ( y k ):˙ ϕ k = ω k + ε ∂ Φ k ∂ y k G k ( y , y , . . . ) . (2)This equation is exact, but it contains the full state space trajectories { y j } . However,for a small perturbation these trajectories are close to the limit cycle, up to deviationsof order ∼ ε . Thus, substituting the zero-order approximation y j ≈ Y j into the r.h.s.of (2), we obtain equations, where only the states on the limit cycles appear, and thesestates are unambiguously determined by the phases:˙ ϕ k = ω k + ε ∂ Φ k ∂ Y k G k ( Y , Y , . . . ) = ω k + εG k ( ϕ k , ϕ , ϕ , . . . ) . As is clear from the presented consideration, to extend the reduction beyond thefirst-order approximation, one has to calculate the deviations from the limit cycle oforders ∼ ε, ε , . . . and to express these deviations in terms of the phases. Below, we donot provide a general approach valid for arbitrary systems, but restrict ourselves to thesimplest case where the phase can be introduced analytically. It is instructive to introduce the Stuart-Landau oscillators in dimensional variables, tounderstand the meaning of dimensionless parameters in the problem. We write theequation for the complex amplitude of oscillations Z for a particular oscillator as dZdτ = ( µ + iν ) Z − ( β + iγ ) | Z | Z + (cid:15)G ( Z , Z , . . . ) , (3)where G is the coupling term depending on the states of other elements of the network.Notice that for brevity and without loss of generality we omit the index for the consideredoscillator and label other units by the index k = 1 , , . . . .In Eq. (3), all the parameters are dimensional, and it is convenient to performa transformation to dimensionless variables and parameters. Here one should alsotake into account, that generally the parameters might be different for differentoscillators. Using only local parameters, one can introduce the dimensionless amplitude Z = (cid:112) µ/βA , so that all uncoupled oscillators will have amplitude one. Another variablethat we can scale is the time, and here one has to choose some common scaling for alloscillators. It appears convenient to assume that the growth rate of linear oscillations µ is the same for all oscillators (the results can be straightforwardly generalized also tothe case of different µ ), and use it to introduce dimensionless time as t = µτ . Then weobtain ˙ A = (1 + iω ) A − | A | A − iαA ( | A | −
1) + εG ( A , A , . . . ) , (4)where ˙ A now means derivative with respect to t . Here ω = ν/µ − γ/β is the dimensionlessfrequency of the limit cycle oscillation and the limit cycle has unit amplitude | A (0) | = 1.Parameter α = γ/β measures non-isochronicity, as we will see below, it determines the igh-Order Phase Reduction for Coupled Oscillators A ). The dimensionless coupling parameter ε depends on thescaling of the function G . If G contains first powers of the amplitudes Z , then ε = (cid:15)/µ ,where we assume that µ/β ≈ µ k /β k , i.e. all the coupled oscillators Z have similaramplitudes. As we see, all the dimensionless parameters entering the problem, namely ω, α, ε , are the original parameters of the system normalized by the linear growth rateof oscillations µ .For further analysis it is instructive to write the Stuart-Landau equations in polarcoordinates, with the amplitude R and the angle θ , so that A = Re iθ and˙ R = R − R + ε Re (cid:2) e − iθ G (cid:3) , (5)˙ θ = ω − α ( R −
1) + εR − Im (cid:2) e − iθ G (cid:3) . (6)It is straightforward to check that, in the absence of coupling, the quantity ϕ = θ − α ln R grows uniformly in time and, hence, is the true phase. Therefore, we rewrite Eqs. (5,6)in terms of R, ϕ :˙ R = R − R + ε Re (cid:2) e − i ( ϕ + α ln R ) G ( R , ϕ , R , ϕ , . . . ) (cid:3) , (7)˙ ϕ = ω + εR [Im (cid:2) e − i ( ϕ + α ln R ) G ( R , ϕ , R , ϕ , . . . ) (cid:3) − α Re (cid:2) e − i ( ϕ + α ln R ) G ( R , ϕ , R , ϕ , . . . ) (cid:3) . (8)Here we have explicitly used that the coupling terms depend on the amplitudes and thephases of other oscillators. Now, we outline the exploited perturbation procedure, while in the next subsection wewill elaborate on it for a particular example. The idea is to look for a dynamics, in whichthe amplitudes are enslaved by the phases. Namely, we assume that the amplitude R isa function of the phases: R = 1 + εr (1) ( ϕ, ϕ , ϕ , . . . ) + ε r (2) ( ϕ, ϕ , ϕ , . . . ) + . . . , (9)and that the dynamics of the phases is also represented as a power series in ε :˙ ϕ = ω + ε Ψ (1) ( ϕ, ϕ , ϕ , . . . ) + ε Ψ (2) ( ϕ, ϕ , ϕ , . . . ) + . . . . (10)We substitute these expressions, together with the similar expressions for R k , ϕ k , intoEqs. (7,8). Equating the terms in each power of ε , we obtain the unknown functions r (1) , r (2) , Ψ (1) , Ψ (2) , . . . .We illustrate here the first few steps. In the equation for the phase, the first-orderapproximation, as discussed above, corresponds to taking into account only the leadingterm in Eq. (9), this yieldsΨ (1) = Im (cid:2) e − iϕ G (1 , ϕ , , ϕ , . . . ) (cid:3) − α Re (cid:2) e − iϕ G (1 , ϕ , , ϕ , . . . ) (cid:3) . The equation for the amplitude in the first order is obtained by substituting Eq. (9) inEq. (7): dr (1) dt = − r (1) + Re (cid:2) e − iϕ G (1 , ϕ , , ϕ , . . . ) (cid:3) . igh-Order Phase Reduction for Coupled Oscillators r (1) via partial derivatives with respect to the phases,and insert the zero-order expressions for these derivatives: dr (1) dt = ∂r (1) ∂ϕ ˙ ϕ + ∂r (1) ∂ϕ ˙ ϕ + ∂r (1) ∂ϕ ˙ ϕ + . . . ≈≈ ∂r (1) ∂ϕ ω + ∂r (1) ∂ϕ ω + ∂r (1) ∂ϕ ω + . . . . (11)Thus, the problem of determining first-order correction to the amplitude reduces to thepartial differential equation2 r (1) + ∂r (1) ∂ϕ ω + ∂r (1) ∂ϕ ω + ∂r (1) ∂ϕ ω + . . . =Re (cid:2) e − iϕ G (1 , ϕ , , ϕ , . . . ) (cid:3) = (cid:88) m,m ,m ,... g mm m ... e i ( mϕ + m ϕ + m ϕ + ... ) (12)Here, we used that the function G is 2 π -periodic with respect to the phases and, hence,the r.h.s. can be written as a multiple Fourier series. Similarly, we expand the function r (1) as: r (1) = (cid:88) m,m ,m ,... ρ mm m ... e i ( mϕ + m ϕ + m ϕ + ... ) , which finally yields an expression for the Fourier coefficients of r (1) ρ mm m ... = g mm m ... imω + im ω + im ω + . . . . (13)These are the basic steps in the perturbation expansion. The expressions for r (1) , r (1) k ,being substituted in Eq. (8) provide phase equations in order ∼ ε . Equations for r (2) are partial differential equations of type (12) with r.h.s containing also Ψ (1) , r (1) , etc. In this Section, we exemplify the perturbative procedure for the derivation of phasedynamics equations in a higher order of the perturbation parameter ε with a particularconfiguration of three coupled Stuart-Landau oscillators. We consider an array of three elements with thecoupling structure 1 ↔ ↔ A = (1 + i ω ) A − | A | A − i αA ( | A | −
1) + εc , e i β , A , ˙ A = (1 + i ω ) A − | A | A − i αA ( | A | −
1) + ε (cid:0) c , e i β , A + c , e i β , A (cid:1) , ˙ A = (1 + i ω ) A − | A | A − i αA ( | A | −
1) + εc , e i β , A . (14)Notice that the oscillators have different frequencies. Introducing the amplitudes andthe angle variables according to A k = R k exp[i θ k ] and the phases ϕ k = θ k − α ln R k , we igh-Order Phase Reduction for Coupled Oscillators R = R − R + εc , R cos( θ − θ + β , ) , ˙ R = R − R + εc , R cos( θ − θ + β , )+ εc , R cos( θ − θ + β , ) , ˙ R = R − R + εc , R cos( θ − θ + β , ) , ˙ ϕ = ω + εc , R R [sin( θ − θ + β , ) − α cos( θ − θ + β , )] , ˙ ϕ = ω + εc , R R [sin( θ − θ + β , ) − α cos( θ − θ + β , )]+ εc , R R [sin( θ − θ + β , ) − α cos( θ − θ + β , )] , ˙ ϕ = ω + εc , R R [sin( θ − θ + β , ) − α cos( θ − θ + β , )] . (15) Now we expand the amplitude deviations in powersof ε (below k = 1 , , R k = 1 + εr (1) k + ε r (2) k + ε r (3) k + . . . . (16)According to this, the angles θ k can be represented as θ k = ϕ k + α (cid:20) εr (1) k + ε r (2) k − . ε (cid:16) r (1) k (cid:17) (cid:21) + . . . . (17)Also the ratios of the amplitudes, entering (15), can be expressed as R m R k = 1 + ε (cid:104) r (1) m − r (1) k (cid:105) + ε (cid:20) r (2) m − r (2) k − r (1) m r (1) k + (cid:16) r (1) k (cid:17) (cid:21) + . . . . (18)Substituting these expansions in Eq. (15), we obtain the following expressions for the igh-Order Phase Reduction for Coupled Oscillators ε :˙ ϕ = ω + εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ ε c , (1 + α ) sin( ϕ − ϕ + β , )( r (1)2 − r (1)1 )+ ε c , (1 + α ) (cid:104) (cid:16) r (2)2 − r (2)1 − r (1)2 r (1)1 + ( r (1)1 ) (cid:17) sin( ϕ − ϕ + β , )+ α (cid:16) . r (1)1 ) + 0 . r (1)2 ) − r (1)2 r (1)1 (cid:17) cos( ϕ − ϕ + β , ) (cid:105) + . . . , ˙ ϕ = ω + εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ ε c , (1 + α ) sin( ϕ − ϕ + β , )( r (1)1 − r (1)2 )+ ε c , (1 + α ) sin( ϕ − ϕ + β , )( r (1)3 − r (1)2 )+ ε c , (1 + α ) (cid:104) (cid:16) r (2)1 − r (2)2 − r (1)1 r (1)2 + ( r (1)2 ) (cid:17) sin( ϕ − ϕ + β , )+ α (cid:16) . r (1)2 ) + 0 . r (1)1 ) − r (1)1 r (1)2 (cid:17) cos( ϕ − ϕ + β , ) (cid:105) + ε c , (1 + α ) (cid:104) (cid:16) r (2)3 − r (2)2 − r (1)3 r (1)2 + ( r (1)2 ) (cid:17) sin( ϕ − ϕ + β , )+ α (cid:16) . r (1)2 ) + 0 . r (1)3 ) − r (1)3 r (1)2 (cid:17) cos( ϕ − ϕ + β , ) (cid:105) + . . . , ˙ ϕ = ω + εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ ε c , (1 + α ) sin( ϕ − ϕ + β , )( r (1)2 − r (1)3 )+ ε c , (1 + α ) (cid:104) (cid:16) r (2)2 − r (2)3 − r (1)2 r (1)3 + ( r (1)3 ) (cid:17) sin( ϕ − ϕ + β , )+ α (cid:16) . r (1)3 ) + 0 . r (1)2 ) − r (1)2 r (1)3 (cid:17) cos( ϕ − ϕ + β , ) (cid:105) + . . . . (19)Next, we have to evaluate corrections r (1) k , r (2) k , . . . . This is accomplished bysubstituting expressions (16) in the equations for the amplitudes in (15). Here the timederivatives are calculated according to the chain rule, as the corrections r (1) k , r (2) k , . . . areassumed to be functions of the phases ϕ k . For the sake of brevity, we present only the igh-Order Phase Reduction for Coupled Oscillators r : ω ∂r (1)1 ∂ϕ + ω ∂r (1)1 ∂ϕ + ω ∂r (1)1 ∂ϕ + 2 r (1)1 = c , cos( ϕ − ϕ + β , ) ,ω ∂r (2)1 ∂ϕ + ω ∂r (2)1 ∂ϕ + ω ∂r (2)1 ∂ϕ + 2 r (2)1 = − r (1)1 ) − αc , ( r (1)2 − r (1)1 ) sin( ϕ − ϕ + β , ) + c , r (1)2 cos( ϕ − ϕ + β , )+ c , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )] ∂r (1)1 ∂ϕ + (cid:104) c , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ c , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )] (cid:105) ∂r (1)1 ∂ϕ + c , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )] ∂r (1)1 ∂ϕ . (20)The partial differential equations for r (1) k , r (2) k , . . . are straightforwardly solved in theFourier representation, because the variations of the amplitudes are 2 π -periodicfunctions of the phases, as outlined in discussion of Eq. (13) above. As a result, weobtain in the 1st order in ε : r (1)1 = 2 c , ω − ω ) cos( ϕ − ϕ + β , ) + ( ω − ω ) c , ω − ω ) sin( ϕ − ϕ + β , ) ,r (1)2 = 2 c , ω − ω ) cos( ϕ − ϕ + β , ) + ( ω − ω ) c , ω − ω ) sin( ϕ − ϕ + β , )+ 2 c , ω − ω ) cos( ϕ − ϕ + β , ) + ( ω − ω ) c , ω − ω ) sin( ϕ − ϕ + β , ) ,r (1)3 = 2 c , ω − ω ) cos( ϕ − ϕ + β , ) + ( ω − ω ) c , ω − ω ) sin( ϕ − ϕ + β , ) . (21)Substitution of these expression in Eq. (19) completes the second-order phasereduction and yields closed equations:˙ ϕ = ω + εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ ε (cid:104) a (2)1;0 + a (2)1; − , , cos(2 ϕ − ϕ ) + b (2)1; − , , sin(2 ϕ − ϕ )+ a (2)1; − , , − cos(2 ϕ − ϕ − ϕ ) + b (2)1; − , , − sin(2 ϕ − ϕ − ϕ )+ a (2)1; − , , cos( ϕ − ϕ ) + b (2)1; − , , sin( ϕ − ϕ ) (cid:105) , (22) igh-Order Phase Reduction for Coupled Oscillators ϕ = ω + εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ ε (cid:104) a (2)2;0 + a (2)2;2 , − , cos(2 ϕ − ϕ ) + b (2)2;2 , − , sin(2 ϕ − ϕ )+ a (2)2;0 , − , cos(2 ϕ − ϕ ) + b (2)2;0 , − , sin(2 ϕ − ϕ )+ a (2)2; − , , − cos(2 ϕ − ϕ − ϕ ) + b (2)2; − , , − sin(2 ϕ − ϕ − ϕ )+ a (2)2;1 , , − cos( ϕ − ϕ ) + b (2)2;1 , , − sin( ϕ − ϕ ) (cid:105) , (23)˙ ϕ = ω + εc , [sin( ϕ − ϕ + β , ) − α cos( ϕ − ϕ + β , )]+ ε (cid:104) a (2)3;0 + a (2)3;0 , , − cos(2 ϕ − ϕ ) + b (2)3;0 , , − sin(2 ϕ − ϕ )+ a (2)3; − , , − cos(2 ϕ − ϕ − ϕ ) + b (2)3; − , , − sin(2 ϕ − ϕ − ϕ )+ a (2)3;1 , , − cos( ϕ − ϕ ) + b (2)3;1 , , − sin( ϕ − ϕ ) (cid:105) . (24)The coefficients of the second-order coupling terms, denoted in Eqs. (22-24) by a (2) k ; l , b (2) k ; l , are listed in Tables A1, A2, A3). Here the 3-component vector l = ( l , l , l )is used to signify the term with the combination of the phases l ϕ + l ϕ + l ϕ .Furthermore, in Table B1 we list the terms (without coupling coefficients) appearing inorders ε , ε .We now shortly discuss the physical meaning of the terms that appear in higherorders in ε .(i) In the second-order approximation, there are no correction terms to the first-ordercouplings. These corrections appear in the third order (and, presumably, in all oddorders).(ii) There are terms, which can be roughly described as “squares” of the basic couplingterms; for the dynamics of the 1st oscillator ϕ these are constant terms and termscontaining the second harmonics of the phase difference, e.g. ∼ sin(2 ϕ − ϕ ).These high-order terms do not arise from the interaction within the whole network;they appear already in a system of two coupled oscillators.(iii) Terms containing combinations of all three phases, e.g., ∼ sin(2 ϕ − ϕ − ϕ ), meanthat an effective hypernetwork with non-pairwise coupling appears already in thesecond-order reduction (cf. studies of synchronization on hypernetworks of phaseoscillators [41, 42]).(iv) Terms containing phase differences for not directly coupled oscillators, (e.g., theterm ∼ sin( ϕ − ϕ ) on the r.h.s. of the equation for ˙ ϕ ) mean that connections interms of phase dynamics do not coincide with the “structural” connections in theoriginal formulation.(v) While the first-order coupling terms are frequency-independent, the second-orderterms depend explicitly on frequency differences. In a more general setup (cf. asystem of coupled van der Pol equations treated below) we expect that couplingwill depend on the frequencies themselves. igh-Order Phase Reduction for Coupled Oscillators
3. Numerical phase reduction
Stuart-Landau oscillator represents an exceptional case when the phase and theinstantaneous frequency can be directly obtained from equations. We exploited thisfeature to derive the second-order phase dynamics equation in the previous Section.For a general oscillator, one has to evaluate the phases numerically. This immediatelyprovides the first-order approximation of the phase dynamics via numerically calculatedphase sensitivity functions.In this Section, we describe a numerical procedure for determining couplingfunctions in higher orders, by virtue of the phase analysis of numerically obtainedtrajectories of the full system. We will first verify it by comparing the results for theStuart-Landau model (14) with theory in Section 2, and then apply it to a network ofthree interacting van der Pol oscillators.
The first step in numerical analysis is the determination of phases ϕ i and their derivatives˙ ϕ i for all elements of the analyzed network, i = 1 , . . . , N . For this purpose, we extendthe technique suggested in [43], where it was described how to obtain phases from anarbitrary trajectory.Consider a particular oscillator within a network described by Eq. (1). Omittingthe index of this unit for simplicity of presentation, we write d y dt = f ( y ) + ε G , (25)where the last term describes the coupling to all other units. As a preparatory step wecompute the autonomous period T . It is, for ε = 0 we take an arbitrary point on thelimit cycle Y , and assign to it ϕ = 0; the return time to this point is exactly T . Now,for any other point on the limit cycle we can compute the time τ ( Y ) required to reachthe zero point where ϕ = 0 or, equivalently, ϕ = 2 π . Since the true phase grows linearlyin time and gains 2 π with one revolution around the cycle, its value can be obtained as ϕ ( Y ) = 2 π T − τ ( Y ) T .The next step is to compute ϕ ( t ) and ˙ ϕ ( t ) for an arbitrary trajectory, on which atsome time t the oscillator has the state u = y ( t ), and the time derivative of this stateis v = ˙ y ( t ). To this end, we introduce an autonomous copy of the investigated unit: d w dt = f ( w ) , (26)and let this auxiliary system evolve from initial conditions w (0) = u , for the time interval nT , where the integer n shall be large enough to ensure that the trajectory attracts tothe limit cycle. (Practically, we stop the evolution when (cid:107) w (( n − T ) − w ( nT ) (cid:107) issmaller than a given tolerance.) Since the time interval of the evolution is a multiple ofthe period, the initial point w (0) = u and the end point w ( nT ) = ¯ w have same valueof the phase. The point ¯ w is on the limit cycle, and, hence, its phase ϕ can be easily igh-Order Phase Reduction for Coupled Oscillators ϕ ( u ) = ϕ ( w (0)) = ϕ ( ¯ w ) = 2 π T − τ ( ¯ w ) T is exactly thedesired phase of the oscillator at the state u .In order to obtain the phase derivative we also have to follow the autonomousevolution (26) towards the limit cycle of the initial condition u + v dt . (Practically,it can be performed simultaneously with the evolution of the point u .) Since dt is(infinitely) small, this can be done by tracing the linear evolution of v dt to ¯ v dt withintime interval nT . The law of this linear evolution is given by the Jacobian of the originalequations (26). Thus, two states u and u + v dt of the coupled system (25) map to twopoints ¯ w and ¯ w + ¯ v dt on the limit cycle of the autonomous system (26). These pointsare characterized by the phases ϕ and ϕ + dϕ , respectively. On the other hand, letus consider the evolution of the autonomous system, from the point ¯ w to ¯ w + ¯ v dt .This evolution occurs along the limit cycle of the system (26) within time interval dt .(Notice that generally dt (cid:54) = dt ). The evolution is governed by the flow on the cycle,i.e. ¯ w + ¯ v dt = ¯ w + f ( ¯ w ) dt , what yields dt = dt (¯ v · f ( ¯ w )) / (cid:107) f ( ¯ w ) (cid:107) . Accordingly, phasegrowth is determined by the natural frequency ω = 2 π/T , i.e. dϕ = ωdt , which finallyyields dϕdt = ω ¯ v · f ( ¯ w ) (cid:107) f ( ¯ w ) (cid:107) . Thus, with the presented algorithm we can compute phases and their derivativesas time series with an arbitrary time step and of sufficient length. Certainly, this canbe done for all elements of the network.
Now, we discuss how the phase dynamics equations of a network can be constructednumerically from given phases and their derivatives. To explain the procedure, it isconvenient to denote the r.h.s. of Eq. (10) as Q k ( ϕ , ϕ , . . . , ϕ N ), where the subscript k = 1 , , . . . , N labels the oscillator. Q k are commonly referred to as the couplingfunctions. Since each Q k is 2 π -periodic with respect to all its arguments, we can writeit as a multiple Fourier series˙ ϕ k = a k ; + (cid:88) l (cid:54) = [ a k ; l cos( ϕ · l ) + b k ; l sin( ϕ · l )] , (27)where l = ( l , l , . . . , l N ) and ϕ = ( ϕ , ϕ , . . . , ϕ N ) are N -dimensional vectors of integerindices and phases, respectively. ϕ · l = (cid:80) Nj =1 l j ϕ j denotes the scalar product betweenthe vectors of phases and the mode indices. Notice the relation between the Fouriercoefficients a k ; l , b k ; l , and coefficients in Eqs. (22-24): a k ; l = a (0) k ; l + εa (1) k ; l + ε a (2) k ; l . . . , b k ; l = εb (1) k ; l + ε b (2) k ; l . . . . Thus, we have time series of the phases ϕ ,n , ϕ ,n , . . . , ϕ N,n and their derivatives˙ ϕ ,n , ˙ ϕ ,n , . . . , ˙ ϕ N,n , where n = 1 , , . . . , L is the point index, Eq. (27) becomes asystem of L linear equations for the unknown coupling coefficients a, b . It is naturalto approximate Q k by a finite Fourier series, preserving only M (cid:28) L terms in the igh-Order Phase Reduction for Coupled Oscillators | l j | ≤ m , whichleaves M = 2( m + 1)(2 m (2 m + 1) + 1) − A k ; l , B k ; l . Thus, as a result of the numerical evaluation we obtainan approximation of the phase dynamics as˙ ϕ k = A k ; + (cid:88) l (cid:54) = , | l j |≤ m [ A k ; l cos( ϕ · l ) + B k ; l sin( ϕ · l )] . (28)The sum in Eq. (27) goes over the combination of mode indices such that either l or − l is counted.The success of this numerical approach depends on how how interdependent theseries φ k are. Here we use two different protocols. Asynchronous case.
Suppose that the network does not synchronize. It means that thetrajectories of the system (1) are dense on the N -dimensional torus. Then one just hasto calculate one long trajectory of the system (1) starting from some initial conditionsand compute phases and instantaneous frequencies for each point of the solution, asdescribed above. If the data series is sufficiently long, the computed trajectory covers thesurface of the torus spanned by ϕ, ϕ , ϕ , . . . , ϕ N . Hence, we obtain enough informationto recover the function Q of N variables and to determine the Fourier coefficients. Synchronous case.
If the network synchronizes, the trajectory on the torus collapsesto a closed line and Eqs. (27) for the Fourier coefficients become dependent and cannotbe solved. However, there is a way to obtain enough data to solve the problem evenin this case. For this goal one considers not one long trajectory, but an ensembleof short transients (cf. [45, 46]). Namely, one starts numerical integration with someasynchronous initial conditions and follows the trajectory unless it is on the invarianttorus. The corresponding transient time should be larger than the amplitude relaxationtime, but smaller than the synchronization time of the phases. Then the procedureis repeated with new initial conditions. Thus, instead of one long record one collectsmany dynamical states on the torus, unless the sufficient number of points is obtained.Certainly, this protocol can be use for the asynchronous case as well.We exemplify these two protocols in the next Section, but before proceeding withthe examples, we have to discuss an important issue. As we mentioned, the least squaresoptimization requires that data points fill the surface of an N -dimensional torus. Thus,we face the curse of dimensionality: the data requirement grows fast with the networksize N and becomes hardly feasible already for N >
3. Some information about thenetwork, e.g., strength of directed links, can, however, be revealed for
N > igh-Order Phase Reduction for Coupled Oscillators Here, we perform the numerical analysis for the system (14) with the goal to verify mainresults in Sec. 2.4 as well as the numerical approach. We choose two sets of parametersthat correspond to asynchronous and synchronous dynamics, respectively, and employthe corresponding protocols.The parameter values are: α = 0 .
1, while ω , ω and ω are varied to changethe synchronization behaviour. Uncoupled oscillators with these parameters have alimit cycle with R = 1 and initial conditions were chosen on it so that the relaxationtime is significantly reduced. The coupling parameters are all equal c j,k = 1 for( j, k ) ∈ { (1 , , (2 , , (2 , , (3 , } and the phase lags are β , = 0 . β , = 0 . β , = 0 .
43 and β , = 0 .
18. We consider two sets of frequencies: case I with ω = −√ / ω = ( √ − /
10 and ω = 0 .
8; here the dynamics is asynchronousquasiperiodic; while for case II with ω = − . ω = 0 and ω = 0 .
33 the dynamicsis synchronous with a long transient. In both cases, we generated the set of points asfollows: we started the dynamics of system (14) with random phases and amplitudesequal to one. Then, after an initial transient time ∆ t = 20 (which has been chosen toensure that the relaxation of the amplitude to the invariant torus is over, but lockingof the phases still does not occur), the values of the phases and their velocities werestored. Altogether we constructed a set of L = 10 data points.Next, we computed the coefficients of the truncated Fourier series, see Eq. (28),for different values of ε , using the SVD approach as described above. The system ofcoupled Stuart-Landau oscillators is invariant with respect to a phase shift ϕ k → ϕ k + φ which means that only modes that fulfill the condition (cid:88) j =1 l j = 0 (29)can exist. To incorporate this into the analysis, we exploited two approaches:(i) only modes satisfying (29) are taken into account, all other modes are set tozero; these modes constitute a small subset of all possible modes and thereforewe determined the Fourier coefficients up to harmonics | l j | ≤ m = 8.(ii) all modes with | l j | ≤ m = 4 are determined.Here we present the results for the case (i), while the results for (ii) are presented inAppendix C.First of all, we compare the theoretical findings presented in Eqs. (22-24) withnumerical results. To this end, we show in Figs. 1,2 the differences | a k ; l − A k ; l | , | b k ; l − B k ; l | , where k = 1 , ,
3, for the theoretically known modes (see Eqs. (22-24)),for the asynchronous and synchronous configurations, respectively.We see that for weak coupling the difference is on the level of numerical precision; itbecomes of the order of one only for such strong coupling as ε = 0 .
4. Next, we see thatthe difference for the first-order terms grows proportionally to ε , in correspondence with igh-Order Phase Reduction for Coupled Oscillators -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 (b) e rr o r( c oup l c oe f ) epsilon10 -10 -8 -6 -4 -2 (a) e rr o r( c oup l c oe f ) ε | a k ; l − A k ; l | , | b k ; l − B k ; l || a k ; l − A k ; l | , | b k ; l − B k ; l | Figure 1.
Results for the asynchronous configuration of three coupled Stuart-Landauoscillators, see Eq. (15). Differences between the numerically calculated couplingcoefficients and the theoretical values given by Eqs. (22-24) are shown as a function ofthe coupling strength, for all oscillators. Panel (a) presents this difference for the termsappearing in the 1st order in ε ; the black dashed line here corresponds to ∼ ε . Panel(b) presents the terms appearing in the 2nd order in ε ; the magenta dashed line herecorresponds to ∼ ε . Red squares (green circles) represent cosine (sine) coefficients.Additional black markers in (b) show the differences between the zero-order terms A k ;0 , k = 1 , , our theoretical conclusion that there are no second-order correction to the first-orderterms. Thus, numerical results exhibit a good correspondence with the theory.Figures 3,4 present an overview of all Fourier coefficients for which we do nothave theoretical values. Namely, we show here all coefficients except for those enteringEqs. (22-24) and analyzed in Figs. 1,2. The results are in full agreement with the power-series representation of the coupling terms. Indeed, we see four groups of coefficientsthat scale as ε , ε , ε and ε , respectively.Finally, we mention that the overall precision of the numerical procedure can beestimated by computing the rest term ξ k of the Fourier series representation in Eq. (28).The results presented in Fig. 5 show that the rest term grows with the coupling strengthas ε . This is an indication that all the terms of orders from 0 to 8 are within the set ofchosen Fourier modes. Up to this point, our analysis was restricted to the case of Stuart-Landau oscillatorswhere expressions for phases and their derivatives are known. In this Section we present igh-Order Phase Reduction for Coupled Oscillators -16 -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 (b) e rr o r( c oup l c oe f ) epsilon10 -10 -8 -6 -4 -2 (a) e rr o r( c oup l c oe f ) ε | a k ; l − A k ; l | , | b k ; l − B k ; l || a k ; l − A k ; l | , | b k ; l − B k ; l | Figure 2.
The same as in Fig. 1, but for the synchronous configuration. ε -18 -16 -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 c oup l c oe f epsilon H k ; l Figure 3.
Results for the asynchronous configuration of three coupled Stuart-Landauoscillators, see Eq. (15). Here amplitudes H k ; l = (cid:113) A k ; l + B k ; l of all Fourier coefficientsfor all three oscillators except for those shown in Fig. 1 are plotted vs. the couplingstrength. Red triangles up, green squares, blue circles, and brown triangles down showthe coefficients that scale as ε , ε , ε , and ε , respectively. (The dashed lines, fromtop to bottom, have slopes 3, 4, 5, and 6 in log-log coordinates.) We did not checkedfor scaling ∼ ε and higher, and show all the coupling coefficients that do not fulfillabove scaling laws with gray diamonds. igh-Order Phase Reduction for Coupled Oscillators ε -18 -16 -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 c oup l c oe f epsilon H k ; l Figure 4.
The same as in Fig. 3, but for the synchronous configuration. -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 e rr o r(r e c on s t r u c t i on ) epsilon ε ξ k Figure 5.
Accuracy of the phase dynamics reconstruction for all oscillators. Redsquares and blue circles represent the rest term ξ k of the Fourier representation Eq. (28)for asynchronous and synchronous cases, respectively. The dashed line shows powerlaw ∼ ε . a purely numerical evaluation of high-order coupling terms for a network of three non-identical van der Pol oscillators:¨ x − µ (1 − x ) ˙ x + ω x = εx ¨ x − µ (1 − x ) ˙ x + ω x = ε ( x + x )¨ x − µ (1 − x ) ˙ x + ω x = εx (30)We fixed parameters ω = 1, ω = 1 . ω = ω (here ω is the spiral mean,the root of cubic equation ω − ω − µ = 1 and varied the couplingconstant ε in the range 0 . ≤ ε ≤ .
3, in this range the dynamics of the network isasynchronous. From a trajectory of the system (30) we obtained time series ϕ , , , ˙ ϕ , , of length L = 10 points each. We used this set to estimate the coefficients of thetruncated Fourier series in Eq. (28) for m = 4. Thus, 729 coupling constants A k ; l , B k ; l were calculated for each ε . Below we restrict our attention only to the strengths of thecoupling and ignore the relative phase, so we look on the properties of 364 couplingconstants H k, l = ( A k ; l + B k ; l ) / . Together with the free term A k ; this constitutes aset of 365 numerically obtained coefficients for each oscillator and for each coupling igh-Order Phase Reduction for Coupled Oscillators ε .For presentation of the results we perform a preliminary sorting. According to thetheory, we expect to obtain terms with leading dependencies ∼ ε q , with q = 0 , , , . . . .Therefore, for each set of indices l and each oscillator, we tried to approximate thecoefficients H k, l ( ε ) by the function ∼ ε q , and if the fitted value q was close to aninteger and the reliability of the fit was high, we attributed the corresponding powerto the coupling term. Additionally, for ε , ε, ε we used the five small values of ε = 0 . , . , . , . , .
02, while for powers ε and ε we used larger couplingconstants ε = 0 . , . , . , . , .
15. We did not look for powers 5 and larger.Figure 6 illustrates these findings. It is instructive to look which coupling modesappear in each power of ε . For modes up to ∼ ε we summarize this in Table 1. Themodes A k ; , describing corrections to the natural frequency are not listed. Notice thatin each cell of the table the modes are ordered according to their amplitudes. Below wesummarize the properties of the coupling modes of the van der Pol oscillators:(i) Looking at the modes that appear in the first order in ε , we notice that theterms with phase differences and the terms with phase sums have nearly the sameamplitude. This means that the coupling terms ∼ ε have nearly the Winfree form:they are products of two functions of the oscillator phase and of the driving oscillatorphase, in full agreement with the first-order theory. Recall that it is different forthe Stuart-Landau oscillators, where only terms with phase differences appear dueto the condition (29).(ii) Because the variable x of the van der Pol equation possesses odd harmonics of thephase (and the same is true for the phase response curve), terms with the thirdharmonics appear already in the 1st order in ε .(iii) Inspection of terms ∼ ε in Table 1 reveals some terms that are not expected inthe first-order analysis, we show them in italic font. Data in Fig. 6 show thatthe amplitudes of these terms are extremely small, comparable to the errors interms reconstruction. We conclude that these terms are probably spurious and justoccasionally possess a scaling ∼ ε , what is not surprising due to the fact that thetotal number of terms to be found is large.(iv) With the size of the time series explored, we could not reliably detect couplingterms appearing in order ε . However, as the figures show, the terms ∼ igh-Order Phase Reduction for Coupled Oscillators -10 -8 -6 -4 -2 -3 -2 -1 (b) c oup l c oe f (cid:1) -10 -8 -6 -4 -2 (a) c oup l c oe f ε H k ; l H k ; l Figure 6.
Coupling coefficients H k, l = ( A k ; l + B k ; l ) / for all three oscillators areshown by red circles, green crosses, and blue pluses, respectively, vs coupling strength ε . Panel (a) shows powers 0 and 1, and panel (b) shows power 2. Dashed black linescorrespond to the scaling ∼ ε and ∼ ε , respectively. -12 -10 -8 -6 -4 -3 -2 -1 (c) c oup l c oe f (cid:1) -12 -10 -8 -6 -4 -2 (b) c oup l c oe f -12 -10 -8 -6 -4 -2 (a) c oup l c oe f ε H k ; l H k ; l H k ; l Figure 7.
The same as Fig. 6, but for powers 3 (a) and 4 (b). Panel (c) presents allother coefficients. Dashed black lines show powers 3,4, and 5, respectively. igh-Order Phase Reduction for Coupled Oscillators ε ε (1,2,0), (3,4,0),(2,1,4)
56 Modes: (0,2,0), (2,-2,0), (1,-2,1), (0,2,-2), (2,0,0),(1,0,1), (1,0,-1), (1,-2,-1), (0,0,2), (1,2,-1), (1,2,1),(2,2,0), (4,-2,0), (0,2,2), (4,0,0), (4,2,0), (0,2,-4),(3,0,1), (3,0,-1), (0,0,4), (1,-4,1), (1,-2,3), (3,2,-1),(1,2,-3), (1,4,-1), (1,-2,-3), (1,0,-3), (0,4,-2), (0,2,4),(1,-4,-1), (3,-2,1), (2,4,0), (3,-2,-1), (1,-4,3), (0,4,0),(1,4,1), (1,2,3), (1,0,3), (2,-4,0), (3,4,-1), (4,-4,0),(0,4,2), (3,2,1), (3,-4,1), (3,0,-3), (3,-2,-3), (0,4,-4),(4,4,0), (1,4,-3), (3,-2,3), (3,0,3), (1,-4,-3), (0,4,4),(1,4,3), (3,-4,-3), (3,2,3)3 9 Modes: (0,1,1),(0,1,-1), (0,3,1), (0,3,-1), (0,1,3), (0,1,-3),(0,3,-3), (0,3,3), (1,-1,1)
50 Modes: (0,2,-2), (0,0,2), (1,0,-1), (1,0,1), (1,-2,-1), (1,-2,1), (0,2,0), (1,2,1), (1,2,-1), (0,2,2), (0,4,-2),(1,-4,-1) , (1,-4,1), (0,4,0), (1,4,1), (1,4,-1), (0,4,2),(3,-2,-1), (3,-4,-1), (1,0,3), (3,-2,1), (3,-4,1), (1,0,-3),(1,-2,-3), (1,-2,3), (0,2,-4), (1,2,3), (1,-4,3), (1,2,-3),(0,2,4) , (3,2,1), (3,2,-1), (3,0,-1), (0,0,4), (3,0,1),(1,4,-3), (1,-4,-3), (1,4,3), (3,-2,-3), (0,4,-4), (3,4,1),(3,-2,3), (0,4,4), (3,4,-1), (3,-4,3), (3,0,3), (3,2,3),(3,2,-3), (3,0,-3), (3,4,3)
Table 1.
All modes revealed in orders ε and ε for a network of van der Pol oscillators.Italic font in the second column denotes the modes that shall not appear in the first-order approximation.
4. Conclusion
In this work, we have presented an analytic perturbation approach, allowing thederivation of the equations for the phase dynamics for general networks of Stuart-Landau oscillators. We exemplified this framework with calculations for a particularsystem of three units. We demonstrated explicitly that already in the second orderin coupling strength, there exist coupling terms that are not present in the structuralcoupling configuration. This result confirms a general statement that phase connectivitygenerally differs from the structural connectivity of a network. In particular, in higher- igh-Order Phase Reduction for Coupled Oscillators
Acknowledgments
E. Gengel acknowledges financial support from the Friedrich-Ebert-Stiftung.A. Pikovsky thanks Russian Science Foundation (Grant Number 17-12-01534) for sup-port. This paper was developed within the scope of the IRTG 1740 / TRP 2015/50122-0,funded by the DFG / FAPESP.
Appendix A. Second-order coupling coefficients in a network ofStuart-Landau oscillators.
Here we present the coefficients for ˙ ϕ , , in Eqs. (22-24). In Tables A1, A2, A3 we usethe notation C m,k = 1 + α c m,k ω m − ω k ) , D m,k = 1 + α ω m − ω k ) c m,k ω m − ω k ) . igh-Order Phase Reduction for Coupled Oscillators a (2)1;0 c , (cid:16) C , sin( β , + β , ) − D , cos( β , + β , ) − D , (cid:17) a (2)1; − , , c , (cid:16) C , sin( β , − β , ) + D , cos( β , − β , ) − C , sin 2 β , + D , cos 2 β , (cid:17) b (2)1; − , , c , (cid:16) C , cos( β , − β , ) − D , sin( β , − β , ) − C , cos 2 β , − D , sin 2 β , (cid:17) a (2)1; − , , − c , (cid:16) C , sin( β , − β , ) + D , cos( β , − β , ) (cid:17) b (2)1; − , , − c , (cid:16) C , cos( β , − β , ) − D , sin( β , − β , ) (cid:17) a (2)1; − , , c , (cid:16) − D , cos( β , + β , ) + C , sin( β , + β , ) (cid:17) b (2)1; − , , c , (cid:16) D , sin( β , + β , ) + C , cos( β , + β , ) (cid:17) Table A1.
Coupling coefficients of the first Stuart-Landau oscillator. a (2)2;0 (cid:16) C , c , sin( β , + β , ) − D , c , cos( β , + β , ) − D , c , + C , c , sin( β , + β , ) − D , c , − D , c , cos( β , + β , ) (cid:17) a (2)2;2 , − , (cid:16) C , c , sin( β , − β , ) + D , c , cos( β , − β , ) − C , c , sin 2 β , + D , c , cos 2 β , (cid:17) b (2)2;2 , − , (cid:16) C , c , cos( β , − β , ) − D , c , sin( β , − β , ) − C , c , cos 2 β , − D , c , sin 2 β , (cid:17) a (2)2;0 , − , (cid:16) C , c , sin( β , − β , ) + D , c , cos( β , − β , ) − C , c , sin 2 β , + D , c , cos 2 β , (cid:17) b (2)2;0 , − , (cid:16) C , c , cos( β , − β , ) − D , c , sin( β , − β , ) − C , c , cos 2 β , − D , c , sin 2 β , (cid:17) a (2)2; − , , − (cid:16) D , c , cos( β , + β , ) − C , c , sin( β , + β , ) − C , c , sin( β , + β , ) + D , c , cos( β , + β , ) (cid:17) b (2)2; − , , − (cid:16) D , c , sin( β , + β , ) + C , c , cos( β , + β , )+ C , c , cos( β , + β , ) + D , c , sin( β , + β , ) (cid:17) a (2)2;1 , , − (cid:16) − D , c , cos( β , − β , ) − C , c , sin( β , − β , ) − C , c , sin( β , − β , ) − D , c , cos( β , − β , ) (cid:17) b (2)2;1 , , − (cid:16) D , c , sin( β , − β , ) − C , c , cos( β , − β , )+ C , c , cos( β , − β , ) − D , c , sin( β , − β , ) (cid:17) Table A2.
Coupling coefficients of the second Stuart-Landau oscillator. igh-Order Phase Reduction for Coupled Oscillators a (2)3;0 c , (cid:16) C , sin( β , + β , ) − D , cos( β , + β , ) − D , (cid:17) a (2)3;0 , , − c , (cid:16) C , sin( β , − β , ) + D , cos( β , − β , ) − C , sin 2 β , + D , cos 2 β , (cid:17) b (2)3;0 , , − c , (cid:16) C , cos( β , − β , ) − D , sin( β , − β , ) − C , cos 2 β , − D , sin 2 β , (cid:17) a (2)3; − , , − c , (cid:16) C , sin( β , − β , ) + D , cos( β , − β , ) (cid:17) b (2)3; − , , − c , (cid:16) C , cos( β , − β , ) − D , sin( β , − β , ) (cid:17) a (2)3;1 , , − c , (cid:16) − D , cos( β , + β , ) + C , sin( β , + β , ) (cid:17) b (2)3;1 , , − c , (cid:16) D , sin( β , + β , ) + C , cos( β , + β , ) (cid:17) Table A3.
Coupling coefficients of the third Stuart-Landau oscillator. igh-Order Phase Reduction for Coupled Oscillators Appendix B. Terms in the higher orders for coupled Stuart-Landauoscillators
In Table B1 we present coupling modes appearing in higher orders in ε . Namely, we justgive the vectors l of these modes. Up to order 4 we checked all of them both analyticallyand numerically; for order 5 we present only numerical results.1 and 3 2 and 4 3 4 51 (-1,1,0) (0,0,0) (-3,3,0), (0,-1,1) (-4,4,0), (0,-2,2) (0,3,-3), (5,-5,0)(-2,2,0) (-1,3,-2), (-2,3,-1) (-2,0,2), (2,-4,2) (2,-5,3), (3,-5,2)(-1,0,1) (-1,-1,2), (-2,1,1) (-1,-2,3), (3,-2,-1) (3,-1,-2), (2,1,-3)(1,-2,1) (-1,4,-3), (-3,4,-1) (4,-5,1), (1,-5,4)(4,-3,-1), (1,3,-4)2 (1,-1,0) (0,0,0) (-3,3,0), (0,3,-3) (-4,4,0), (0,4,-4) (5,-5,0), (0,5,-5)(0,-1,1) (2,-2,0) (-1,3,-2), (-2,3,-1) (-2,0,2), (2,-4,2) (3,-1,-2), (1,2,-3)(0,-2,2) (-1,-1,2), (-2,1,1) (-1,-2,3), (3,-2,-1) (1,3,-4), (4,-3,-1)(-1,0,1) (-1,4,-3), (-3,4,-1) (4,-5,1), (1,-5,4)(1,-2,1) (2,-5,3), (3,-5,2) (3,-5,2), (2,-5,3)3 (0,1,-1) (0,0,0) (0,3,-3), (1,-1,0) (0,4,-4), (2,-2,0) (3,-3,0), (0,5,-5)(0,2,-2) (-1,3,-2), (-2,3,-1) (2,0,-2), (2,-4,2) (3,-5,2), (2,-5,3)(1,0,-1) (-1,-1,2), (-2,1,1) (-1,-2,3), (3,-2,-1) (3,-1,-2), (2,1,-3)(1,-2,1) (-1,4,-3), (-3,4,-1) (1,-5,4), (4,-5,1)(1,3,-4), (4,-3,-1) Table B1.
Coupling terms that appear in different orders (see columns) and for allthree oscillators (see rows).
Appendix C. Numerical reconstruction of coupling for Stuart-Landauoscillators: approach (ii)
Here we present the result of the numerical procedure using the approach (ii), i.e. herewe do not exclude zero modes that do not fulfill the condition (29). Comparison of thiscase with the results for the approach (i), presented in the main text, is important forthe analysis of networks of van der Pol or other oscillators, for which no modes canbe excluded a priori . Figures C1-C5 shall be compared to Figs. 1-5, respectively. Thecomparison shows that the results are consistent. As expected, the approach (i) providesa higher accuracy since fewer unknowns shall be found. However, even with the secondapproach, we managed to reliably reveal scaling of the mode coefficients up to the order ∼ ε . igh-Order Phase Reduction for Coupled Oscillators -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 (b) e rr o r( c oup l c oe f ) epsilon10 -10 -8 -6 -4 -2 (a) e rr o r( c oup l c oe f ) | a k ; l − A k ; l | , | b k ; l − B k ; l || a k ; l − A k ; l | , | b k ; l − B k ; l | ε Figure C1.
Results for the asynchronous configuration (case I) of three coupledStuart-Landau oscillators, see Eq.(15). Differences between the numerically calculatedcoupling coefficients and the theoretical values given by Eqs.(22-24) are shown as afunction of the coupling strength, for all oscillators. Panel (a) presents this differencefor the terms appearing in the 1st order in ε ; the black dashed line here correspondsto ∼ ε . Panel (b) presents the terms appearing in the 2nd order in ε ; the magentadashed line here corresponds to ∼ ε . Red squares (green circles) represent cosine(sine) coefficients. Additional black markers in (b) show the differences between thezero-order terms A k ;0 , k = 1 , , -16 -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 (b) e rr o r( c oup l c oe f ) epsilon10 -10 -8 -6 -4 -2 (a) e rr o r( c oup l c oe f ) ε | a k ; l − A k ; l | , | b k ; l − B k ; l || a k ; l − A k ; l | , | b k ; l − B k ; l | Figure C2.
The same as in Fig. C1, but for the synchronous configuration. igh-Order Phase Reduction for Coupled Oscillators -18 -16 -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 c oup l c oe f epsilon ε H k ; l Figure C3.
Results for the asynchronous configuration (case I) of three coupledStuart-Landau oscillators, see Eq. (15). Here all Fourier coefficients for all threeoscillators except for those shown in Fig. C1 are plotted vs. the coupling strength.Red triangles up, green squares, blue circles, and brown triangles down show thecoefficients that scale as ε , ε , and ε , respectively. (The dashed lines, from top tobottom, have slopes 3, 4, and 5, in log-log coordinates.) We did not checked for scaling ∼ ε and higher, and show all the coupling coefficients that do not fulfill above scalinglaws with gray pluses. -15 -10 -5 -3 -2 -1 c oup l c oe f epsilon ε H k ; l Figure C4.
The same as in Fig. C3, but for the synchronous configuration. -14 -12 -10 -8 -6 -4 -2 -3 -2 -1 e rr o r(r e c on s t r u c t i on ) epsilon ε ξ k Figure C5.
Accuracy of the phase dynamics reconstruction. Red squares and bluecircles illustrate asynchronous and synchronous cases, respectively. The dashed lineshows power law ∼ ε . igh-Order Phase Reduction for Coupled Oscillators [1] Micha Nixon, Eitan Ronen, Asher A. Friesem, and Nir Davidson. Observing geometric frustrationwith thousands of coupled lasers. Phys. Rev. Lett. , 110:184102, 2013.[2] M. H. Matheny, M. Grau, L. G. Villanueva, R. B. Karabalin, M. C. Cross, and M. L. Roukes. Phasesynchronization of two anharmonic nanomechanical oscillators.
Phys. Rev. Lett. , 112:014101,2014.[3] Balth. van der Pol and J. van der Mark. The heartbeat considered as a relaxation oscillation, andan electrical model of the heart.
The London, Edinburgh, and Dublin Philosophical Magazineand Journal of Science , 6(38):763–775, 1928.[4] I Ashraf, R Godoy-Diana, J Halloy, B Collignon, and B Thiria. Synchronization and collectiveswimming patterns in fish (hemigrammus bleheri).
Journal of the Royal Society Interface ,13(123):20160734, 2016.[5] Steven H Strogatz, Daniel M Abrams, Allan McRobie, Bruno Eckhardt, and Edward Ott.Theoretical mechanics: Crowd synchrony on the Millennium Bridge.
Nature , 438(7064):43–44, 2005.[6] Florian D¨orfler and Francesco Bullo. Synchronization in complex networks of phase oscillators: Asurvey.
Automatica , 50(6):1539–1564, 2014.[7] J. Garc´ıa-Ojalvo, J. Casademont, and M. C. Torrent. Coherence and synchronization in diode-laser arrays with delayed global coupling.
Int. J. of Bifurcation and Chaos , 9(11):2225–2229,1999.[8] M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Salova, M. Rohden, J. Li, M. Hudobade Badyn, M. P’osfai, L. Duenas-Osorio, M. Mesbahi, J. P. Crutchfield, M. C. Cross, R. M.DSouza, and M. L. Roukes. Exotic states in a simple network of nanoelectromechanicaloscillators.
Science , 363:eaav7932, 2019.[9] A. B. Cawthorne, P. Barbara, S. V. Shitov, C. J. Lobb, K. Wiesenfeld, and A. Zangwill.Synchronized oscillations in Josephson junction arrays: The role of distributed coupling.
Phys.Rev. B , 60:7575–7578, 1999.[10] V. Tiberkevich, A. Slavin, E. Bankowski, and G. Gerhart. Phase-locking and frustration in anarray of nonlinear spin-torque nano-oscillators.
Appl. Phys. Lett. , 95:262505, 2009.[11] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa. Spontaneous synchrony in power-gridnetworks.
Nat. Phys. , 9:191, 2013.[12] Peter Uhlhaas, Gordon Pipa, Bruss Lima, Lucia Melloni, Sergio Neuenschwander, Danko Nikoli,and Wolf Singer. Neural synchrony in cortical networks: history, concept and current status.
Frontiers in Integrative Neuroscience , 3:17, 2009.[13] L. Glass. Synchronization and rhythmic processes in physiology.
Nature , 410:277–284, 2001.[14] A. Prindle, P. Samayoa, I. Razinkov, T. Danino, L. S. Tsimring, and J. Hasty. A sensing array ofradically coupled genetic “biopixels”.
Nature , 481(7379):39–44, 2012.[15] Yoshiki Kuramoto.
Chemical oscillations, turbulence and waves . Springer, Berlin, 1984.[16] Arkady Pikovsky, Jurgen Kurths, Michael Rosenblum, and J¨urgen Kurths.
Synchronization: auniversal concept in nonlinear sciences . Cambridge University Press, 2001.[17] Hiroya Nakao. Phase reduction approach to synchronisation of nonlinear oscillators.
ContemporaryPhysics , 57(2):188–214, 2016.[18] Bastian Pietras and Andreas Daffertshofer. Network dynamics of coupled oscillators and phasereduction techniques.
Physics Reports , 2019.[19] Y. Kuramoto. Self-entrainment of a population of coupled nonlinear oscillators. In H. Araki,editor,
International Symposium on Mathematical Problems in Theoretical Physics , page 420,New York, 1975. Springer Lecture Notes Phys., v. 39.[20] J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ritort, and R. Spigler. The Kuramoto model:A simple paradigm for synchronization phenomena.
Rev. Mod. Phys. , 77(1):137–175, 2005.[21] Iv´an Le´on and Diego Paz´o. Phase reduction beyond the first order: The case of the mean-fieldcomplex Ginzburg-Landau equation.
Physical Review E , 100(1):012211, 2019.[22] Dan Wilson and Jeff Moehlis. Isostable reduction of periodic orbits.
Physical Review E , igh-Order Phase Reduction for Coupled Oscillators Progress of Theoretical Physics , 88(6):1213–1218, 1992.[24] Michael Rosenblum and Arkady Pikovsky. Self-organized quasiperiodicity in oscillator ensembleswith global nonlinear coupling.
Physical Review Letters , 98(6):064101, 2007.[25] Wataru Kurebayashi, Sho Shirasaka, and Hiroya Nakao. Phase reduction method for stronglyperturbed limit cycle oscillators.
Physical Review Letters , 111(21), November 2013.[26] Kestutis Pyragas and Viktor Noviˇcenko. Phase reduction of a limit cycle oscillator perturbed bya strong amplitude-modulated high-frequency force.
Physical Review E , 92(1), July 2015.[27] Erik Gengel and Arkady Pikovsky. Phase demodulation with iterative Hilbert transformembeddings.
Signal Processing , 165:115–127, 2019.[28] Rok Cestnik and Michael Rosenblum. Inferring the phase response curve from observation of acontinuously perturbed oscillator.
Scientific reports , 8(1):1–10, 2018.[29] Bj¨orn Kralemann, Laura Cimponeriu, Michael Rosenblum, Arkady Pikovsky, and Ralf Mrowka.Phase dynamics of coupled oscillators reconstructed from data.
Physical Review E , 77(6), June2008.[30] B. Bezruchko, V. Ponomarenko, M. G. Rosenblum, and A. S. Pikovsky. Characterizing directionof coupling from experimental observations.
CHAOS , 13(1):179–184, 2003.[31] I. T. Tokuda, S. Jain, I. Z. Kiss, and J. L. Hudson. Inferring phase equations from multivariatetime series.
Phys. Rev. Lett. , 99:064101, 2007.[32] K. A. Blaha, A. Pikovsky, M. Rosenblum, M. T. Clark, C. G. Rusin, and J. L. Hudson.Reconstruction of two-dimensional phase dynamics from experiments on coupled oscillators.
Phys. Rev. E , 84:046201, 2011.[33] B. Kralemann, M. Fr¨uhwirth, A. Pikovsky, M. Rosenblum, T. Kenner, J. Schaefer, and M. Moser.
In vivo cardiac phase response curve elucidates human respiratory heart rate variability.
NatureCommunications , 4:2418, 2013.[34] M. Rosenblum, M. Fr¨uhwirth, M. Moser, and A. Pikovsky. Dynamical disentanglement in ananalysis of oscillatory systems: an application to respiratory sinus arrhythmia.
PhilosophicalTransactions of the Royal Society A: Mathematical, Physical and Engineering Sciences ,377(2160):20190045, 2019.[35] Valentina Ticcinelli, Tomislav Stankovski, Dmytro Iatsenko, Alan Bernjak, Adam E Bradbury,Andrew R Gallagher, Peter Clarkson, Peter VE McClintock, and Aneta Stefanovska. Coherenceand coupling functions reveal microvascular impairment in treated hypertension.
Frontiers inphysiology , 8:749, 2017.[36] C¸ a˘gda¸s Top¸cu, Matthias Fr¨uhwirth, Maximilian Moser, Michael Rosenblum, and Arkady Pikovsky.Disentangling respiratory sinus arrhythmia in heart rate variability records.
Physiologicalmeasurement , 39(5):054002, 2018.[37] Tomislav Stankovski, Tiago Pereira, Peter VE McClintock, and Aneta Stefanovska. Couplingfunctions: universal insights into dynamical interaction mechanisms.
Reviews of ModernPhysics , 89(4):045001, 2017.[38] Tomislav Stankovski, Spase Petkoski, Johan Raeder, Andrew F Smith, Peter VE McClintock, andAneta Stefanovska. Alterations in the coupling functions between cortical and cardio-respiratoryoscillations due to anaesthesia with propofol and sevoflurane.
Philosophical Transactions of theRoyal Society A: Mathematical, Physical and Engineering Sciences , 374(2067):20150186, 2016.[39] Bj¨orn Kralemann, Arkady Pikovsky, and Michael Rosenblum. Reconstructing phase dynamics ofoscillator networks.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 21(2):025104,2011.[40] Bj¨orn Kralemann, Arkady Pikovsky, and Michael Rosenblum. Reconstructing effective phaseconnectivity of oscillator networks from observations.
New Journal of Physics , 16(8):085013,August 2014.[41] M. Komarov and A. Pikovsky. Finite-size-induced transitions to synchrony in oscillator ensembles igh-Order Phase Reduction for Coupled Oscillators with nonlinear global coupling. Phys. Rev. E , 92:020901, 2015.[42] C. C. Gong and A. Pikovsky. Low-dimensional dynamics for higher-order harmonic, globallycoupled phase-oscillator ensembles.
Phys. Rev. E , 100:062210, 2019.[43] Michael Rosenblum and Arkady Pikovsky. Numerical phase reduction beyond the first orderapproximation.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 29(1):011105, 2019.[44] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery.
Numericalrecipes: The art of scientific computing (3rd edition) . Cambridge university press, 2007.[45] Z. Levnaji´c and A. Pikovsky. Network reconstruction from random phase resetting.
Phys. Rev.Lett. , 107(3):034101, 2011.[46] A. Pikovsky. Reconstruction of a random phase dynamics network from observations.
PhysicsLetters A , 382(4):147 – 152, 2018.[47] Bj¨orn Kralemann, Arkady Pikovsky, and Michael Rosenblum. Detecting triplet locking by tripletsynchronization indices.
Physical Review E , 87(5), May 2013.[48] Hannes Osterhage, Florian Mormann, Tobias Wagner, and Klaus Lehnertz. Measuring thedirectionality of coupling: phase versus state space dynamics and application to EEG timeseries.
International journal of neural systems , 17(03):139–148, 2007.[49] Thorsten Rings and Klaus Lehnertz. Distinguishing between direct and indirect directionalcouplings in large oscillator networks: Partial or non-partial phase analyses?