Higher-Order Network Interactions through Phase Reduction for Oscillators with Phase-Dependent Amplitude

Coupled oscillator networks provide mathematical models for interacting periodic processes. If the coupling is weak, phase reduction -- the reduction of the dynamics onto an invariant torus -- captures the emergence of collective dynamical phenomena, such as synchronization. While a first-order approximation of the dynamics on the torus may be appropriate in some situations, higher-order phase reductions become necessary, for example, when the coupling strength increases. However, these are generally hard to compute and thus they have only been derived in special cases: This includes globally coupled Stuart--Landau oscillators, where the limit cycle of the uncoupled nonlinear oscillator is circular as the amplitude is independent of the phase. We go beyond this restriction and derive second-order phase reductions for coupled oscillators for arbitrary networks of coupled nonlinear oscillators with phase-dependent amplitude, a scenario more reminiscent of real-world oscillations. We analyze how the deformation of the limit cycle affects the stability of important dynamical states, such as full synchrony and splay states. By identifying higher-order phase interaction terms with hyperedges of a hypergraph, we obtain natural classes of coupled phase oscillator dynamics on hypergraphs that adequately capture the dynamics of coupled limit cycle oscillators.


Introduction
Collective behavior of oscillatory systems, such as synchronization, is ubiquitous in many realworld dynamical systems.Examples range from biological systems, such as the continuous beat of a healthy heart, regularly spiking neurons, or flashing fireflies to large scale technological systems, such as power grids, or the periodic motion of planets [35,14,22].If the coupling is sufficiently weak, phase reductions provide a useful tool to study these oscillator networks [34,32,33] and have found application to elucidate collective dynamics for example in neuroscience [4].Specifically, if one assumes that each oscillatory unit has a stable limit cycle without coupling, then the whole system still settles to an invariant torus for sufficiently small coupling strengths.A phasereduced system then describes the dynamics on this invariant torus and it is usually derived as an expansion in the coupling strength.Typically, one considers a first-order truncation ignoring any higher-order terms.However, this is insufficient to describe the dynamics of the full/unreduced system when the first-order truncation undergoes a bifurcation or when the coupling is stronger.
To accurately describe the dynamics of the full system in these cases, one has to resort to higher-order phase reductions.Recently progress has been made to compute such phase reductions: Explicit computations show how nonpairwise phase interactions enter the phase-reduced equations once one goes to second or higher orders [31,21].In general, however, computing higher-order phase reductions is not straightforward and the focus has been on simple oscillator models with additive interactions.For example, the computations in [31] make explicit use of the rotational symmetry of Stuart-Landau oscillators that imply that the (unperturbed) limit cycle is the unit circle; the resulting phase equations reflect the symmetry properties.Such an assumption of rotational symmetry is rarely satisfied for general oscillator systems.Indeed, while a limit cycle may be approximately circular for oscillations emanating from a Hopf bifurcation point, oscillations further away from the bifurcation point will generically have a phase-dependent amplitude.
Here, we derive higher-order phase reductions for systems in which the limit cycle has phasedependent amplitude.More specifically, we generalize recent approaches for coupled Stuart-Landau oscillators on a graph (i.e., coupling is additive) to oscillators subject to a perturbation that breaks the rotational symmetry of the limit cycle.We derive phase equations by expanding in terms of both the coupling strength between oscillators as well as the size of the symmetrybreaking perturbation.We then analyze how these higher-order interaction terms affect the stability full synchrony-all oscillators are at the same state-and the splay configuration, in which the phases of the oscillators are uniformly spread out.Our approach allows to compute phase reductions not only for symmetric all-to-all coupled networks, but also coupled oscillators on arbitrary graphs.Thus, for coupled oscillators on a given graph with additive interactions, we obtain a parameterized family of effective phase dynamics which include higher-order nonpairwise phase interactions that depend nonlinearly on three or more oscillator phases.
Our results also provide a tool to construct phase dynamics on hypergraphs that are a meaningful approximation of systems of nonlinearly coupled oscillators.Indeed, nonpairwise phase coupling terms between more than two oscillator phases have been associated with phase oscillator dynamics on hypergraphs [6,9]; these can arise for example through higher-order phase reductions of additively coupled systems (e.g., [31]) or first-order phase reductions of oscillators with generic nonlinear coupling [4].So far, many phase oscillator networks with higher-order interactions that have been considered were ad-hoc, for example, by generalizing the Kuramoto model to hypergraphs (see, e.g., [38,8]).By contrast, our results provide a natural family of hypergraphs together with phase interaction functions that describe the synchronization behavior of an (unreduced) nonlinear oscillator network.This family is parameterized in terms of the underlying coupling graph as well as the system parameters.
This paper is organized as follows: In Section 2 we recall the main points of a previous work [31] considering higher-order phase reductions for globally coupled Stuart-Landau oscillators.Then, in Section 3 we study a system whose limit cycle can be obtained from perturbing the circular limit cycle from a Stuart-Landau oscillator.We derive phase reductions as an expansion in the coupling strength up to second order and the parameter that controls the deformation of the limit cycle.Next, in Section 4 we numerically analyze how the deformation of the limit cycle affects the stability of synchronized and splay states.We investigate how accurately first and second order phase reductions reproduce these stability properties.Finally, Section 6 contains a discussion and some concluding remarks.

Phase Reductions for Stuart-Landau Oscillators
In this section, we recall the main aspects of how to derive higher-order phase reductions for coupled Stuart-Landau oscillators from [31].We highlight the main assumptions that are made to derive these reductions.
First, let us consider a single complex Stuart-Landau oscillator with state A = A(t) ∈ C, that evolves according to where c 2 ∈ R is a parameter.The right-hand side of (2.1) is equivariant with respect to the continuous group S := R/(2πZ), which acts on C by shifting an oscillator by a given phase.Due to this symmetry, it makes sense to introduce polar coordinates A = re iφ with r ≥ 0 and φ ∈ S.Then, the system has a stable limit cycle at r = 1 and the dynamics on this limit cycle can be described by just the phase φ.In fact, on the limit cycle, the phase φ increases with constant speed −c 2 , such that φ(t) = φ(0) − c 2 t.To understand the dynamics off the limit cycle, we note that every point in the basin of attraction of this attractive limit cycle has an asymptotic phase, which is the phase of an initial condition condition of a trajectory on the limit cycle that the point converges to.The set of points with the same asymptotic phase are called isochrons [25,30].To define them, we introduce the notation Φ t A 0 for the solution of (2.1) with A(0) = A 0 at time t.
If A 0 is on the limit cycle, i.e., |A 0 | = 1, and φ(t) = θ − c 2 t then is the isochron with asymptotic phase θ.Upon variation of θ, the isochrons foliate the basin of attraction of the limit cycle.For (2.1), they can explicitly be calculated [31] to be As one can see from this formula these isochrons are symmetric in the sense that one isochron can be obtained from another by shifting it by a constant phase.This property also follows directly from the S symmetry of (2.1).
Having studied a single Stuart-Landau oscillator, we now consider N ∈ N coupled Stuart-Landau oscillators described by complex variables A k , k = 1, . . ., N .When the coupling is as in [31], they satisfy Here, c 1 , c 2 are two real parameters, ǫ ≥ 0 relates to the coupling strength and Ā = 1 N N j=1 A j .Now one changes to polar coordinates A k = r k e iφ k , with r k ≥ 0 and φ k ∈ S := R/(2πZ).After conducting an additional nonlinear transformation θ = φ − c 2 ln(r) to straighten the isochrons, the authors of [31] arrive at the system where ω ∈ R is a parameter that depends on c 2 , f (r) = r(1 − r 2 ) and the functions g k and h k are given by The coupling in (2.2) respects the S symmetry of a Stuart-Landau oscillator such that (2.2) again possesses an S symmetry group that acts on C N by shifting all oscillators A 1 , . . ., A N by the same phase.Since the transformation that straightens the isochrons does not break this symmetry, the system (2.3) inherits the same symmetry group.This fact can also be observed by directly looking at the structure of the functions g k , h k and f .Since they only depend on phase differences, they are invariant when all oscillators are shifted by the same phase.Therefore, we can change into a co-rotating coordinate frame and thereby set ω = 0 without loss of generality.Next, we derive first and second order phase reductions for the system (2.3).In absence of the coupling, i.e., ǫ = 0, each oscillator in (2.2) and (2.3) has a stable limit cycle at |A k | = 1 or r k = 1, respectively.Therefore, the limiting dynamics of the whole system takes place on the N -dimensional torus that is described by r k ≡ 1 for all k = 1, . . ., N .When slightly increasing ǫ this torus persists but it gets perturbed.The radii of this invariant torus are then functions of the phases.In fact, they can be expanded in terms of ǫ such that with r (0) k ≡ 1, see [31].Inserting the ansatz (2.5) into (2.3b) as done in [31] yields θk = ǫh k (r (0) , θ) since ω = 0.By truncating terms of order O(ǫ 2 ) and higher orders, one can obtain a first-order phase reduction.To derive a second order phase reduction one also needs to know r (1) (θ).This can be accomplished by inserting (2.5) into (2.3a) and collecting terms of order ǫ, see [31].One then ends up with Moreover, by the chain rule, one has for the dynamics on the invariant torus, where ½ = (1, . . ., 1) T ∈ R N .Now, it is crucial that ω can be set to 0, because then collecting terms of order ǫ in this equation yields ) and (2.8), the authors of [31] arrive at Substituting that into (2.6) and truncating O(ǫ 3 ) terms yields the second order phase reduction.
3 Phase Reductions for Limit Cycles with Phase-Dependent Amplitude In this section, we introduce a variation of Stuart-Landau oscillators where the limit cycle is not circular but has a phase dependent amplitude.We then derive first-and second-order phase reductions reductions for this class of oscillators subject to coupling as in the previous section.Finally, we investigate how these phase reductions are affected by the parameter that determines the deviation of the shape of the limit cycle from a circle.
Inspired by [31], we start with a modified Stuart-Landau oscillator that can have a noncircular limit cycle of a given functional form.Specifically, for a parameter δ ∈ R with |δ| ≪ 1 and a given smooth 2π-periodic function g : S → R, the oscillator we consider has a limit cycle with phase-dependent amplitude r = 1 + δg(φ).This is the case if the state of oscillator A = re iφ ∈ C evolves according to where ω > 0 is the angular velocity of the oscillator and m < 0 determines the rate of attraction to the limit cycle.The shape of the limit cycle is shown in Figure 1; for δ = 0 the limit cycle is circular.Due to the explicit embedding of the limit cycle, we choose the nonlinearity of the radial direction to be slightly different as in the Stuart-Landau oscillator (2.1).As this primarily influences the rate of radial convergence, this should be not relevant for the phase reductions since the derivatives of the radial equation at r = 1 coincide if m = −2.Moreover, to keep the problem tractable, we focus on oscillations without radial dependency of the phase dynamics (corresponding to c 2 = 0 in (2.1)).Now, given an ensemble of N oscillators (A k = r k e iφ k ) k=1,...,N , we assume a mean-field coupling, of the form where F (A k ) denotes the intrinsic dynamics of oscillator k as described in (3.1),K ∈ R is the coupling strength, α ∈ S is a parameter and Ā = 1 N N j=1 A j is the average position as before.Rewritten in polar coordinates, this results in the system

After the transformation
to transform the phase-dependent limit cycle to a circle, we arrive by a direct, but slightly lengthy, calculation at the system with functions F, G k and H k defined by It is important to note, that for a general perturbation F, G and H depend explicitly on the phase.Therefore, the system (3.3)does not have an S symmetry.Hence, we cannot change into a co-rotating coordinate system.That means, unlike in the system (2.3), we cannot assume ω = 0 without loss of generality.However, ω = 0 was a crucial assumption to derive formula (2.8) in Section 2. Next, we show how to generalize the methods of [31] to derive higher-order phase reductions anyway.

Phase Reductions of First-Order
Because there is an additional parameter δ in our system, we expand in both K and δ, i.e., we study asymptotics in two parameters [29].Regarding notation, if W is a function or a scalar, we write W (n,j) for the contribution of order K n δ j to W , i.e., K n δ j W (n,j) .
Moreover, we write W (n,⋆) for all contributions of order K n , which includes all orders in δ.
Similarly, W (⋆,j) includes all terms of order δ j .Consequently, In particular, if the quantity W is independent of K, we have W = W (0,⋆) , but we use the notation W = W (−,⋆) to highlight the independence of K. We use this notation to derive phase reductions of different approximation order in K and δ.To distinguish these phase reductions, we speak of an (a, b)-phase reduction when a is the highest approximation order in K and b is the highest order in δ.In particular, an (a, b)-phase reduction is given by where P (n,j) k (φ) denotes the contribution on the order K n δ j .Explicit expressions for P (n,j) k (φ) will be derived below.
By the same reasons as illustrated in Section 2, the limiting dynamics of the system (3.3)takes place on an attractive invariant N -dimensional torus.If K = 0 this torus is described by R k ≡ 1 for all k = 1, . . ., N .If |K| is small, the torus persists and the radii of this torus can be expanded in terms of K as where When inserting the expansion (3.4) into the system (3.3b)we obtain which is the base equation for phase reductions of any order.A phase reduction of first order can be obtained by truncating terms of order O(K 2 ), a second order phase reduction is derived from (3.5) by ignoring all terms of order O(K 3 ), etc.In particular, the (1, ∞)-phase reduction is given by Up to now, this contains all orders of δ, but by writing the (1, ∞)-phase reduction can also be written as where (φ) ≡ 0 for all j ∈ N and For example, we find Consequently, the (1, 0)-phase reduction is φk = which is the Kuramoto-Sakaguchi model [36] for identical oscillators.

Higher-Order Phase Reductions
Having explicitly stated the first-order phase reductions, we move on by considering the terms of order K 2 in (3.5), thereby deriving a second-order phase reduction.As in Section 2, this requires knowledge of R (1,⋆) (φ).To get a formula describing it, we follow the lines of [31] and insert the expansion (3.4) into (3.3a).Using R (0,⋆) ≡ 1 and applying the chain rule, we find that the left-hand side of (3.3a) turns into whenever the dynamics is constrained to the limiting torus.Using, Similarly, the right-hand side turns into where ) and (3.8) and collecting terms of order O(K 1 ) yields or equivalently, when using definitions of F and G, which is a linear first-order partial differential equation describing R (1,⋆) k (φ).At this point, we can no longer proceed as in Section 2, because we cannot set ω = 0, since our system is not rotationally invariant.Thus, we generalize the methods of [31] by solving the PDE (3.9), as proposed in [21].Assuming an expansion we solve the PDE (3.9) order by order [18,17,27].When δ = 0, the PDE describing R The solution to this PDE, which can be, for example, be found with the method of characteristics [18], is given by On first-order in δ, the resulting PDE is The solution of this PDE now depends on the specific choice of g.However, as one can infer from the structure of (3.12), its solutions are linear in g in the sense that if R(1,1) is the solution to (3.12) when g = ĝ + g.If g(φ) = sin(φ), the solution of (3.12) is given by where s 1 (φ k , φ l ) is a trigonometric polynomial that is defined by Equivalently, if α = 0, this can also be written as Finally, the PDE on order O(δ 2 ) is This PDE, however, is not linear in g.
where s 2 (φ k , φ l ) is a trigonometric polynomial of the same form as s 1 but with more summands1 .
To determine the second order interactions in (3.5), we also need to expand where e k is the kth unit vector in R N .Finally, we can put everything together and calculate the second order terms in (3.5): Evaluating these expressions yields, for example, for the second-order terms in K when δ = 0. We emphasize that it is here clearly visible that three phases interact with each other, which is different from the terms that appear in a (1, 0) phase reduction.
In conclusion, the (2, 2)-phase reduction is given by φk

Comparison of Phase Reductions with and without Symmetry
As we have highlighted in this section, the full system (3.3) has a rotational S symmetry for δ = 0 that breaks for a generic perturbation to the limit cycle and δ = 0.
We now consider the reduced phase equations in detail focusing on the corresponding change in symmetry properties one would expect2 .The (2, 0)-phase reduction, i.e., the second order phase reduction when δ = 0, is given by As one can see, the right-hand side of this equation only depends on phase differences.Therefore, its value remains invariant when shifting all oscillators by a common phase.Consequently, the (2, 0)-phase reduction inherits the S symmetry of the full system.Writing Qe iΘ = 1 N N j=1 e 2iφj and Re iΨ = 1 N N j=1 e iφj as in [31], we can compare the results in [31] with our (2, 0)-phase reduction.With this notation the (2, 0)-phase reads which agrees with the result from [31, Equation ( 15)], if one chooses This shows that even though the nonlinearity in the radial direction in (3.1a) is different from the nonlinearity of the Stuart-Landau oscillator considered in [31], the phase equations up to second order in K agree as expected.The plausibility of K = ǫ|1 + ic 1 | is immediate, when comparing the coupling strength in (2.2) with those in (3.2).Moreover, m = −2 can be explained as follows: When δ = 0, m is the rate of attraction towards the limit cycle of a single oscillator (3.1).In particular, this rate can be obtained by linearizing (3.1a) with respect to r and evaluating at the limit cycle r = 1.Doing the same for the oscillator (2.1), yields that the rate of attraction to this limit cycle is −2.Now, let us consider phase reductions when δ = 0.When g(φ) = sin(φ), the (2, 1)-phase reduction is given by the system which does not have an S symmetry.In particular, the terms of order δ are not invariant when one shifts all oscillators by a common phase.Since, higher-order phase reductions also consist of these terms, any phase reduction of higher-order is not S symmetric.To conclude, phase reductions of the full system (3.3)possess an S symmetry if and only if the full system itself possesses this symmetry.
As a remark, when α = 0, the (2, 1)-phase reduction can also be written as where This formula shows that effect of δ on the (2, 1)-phase reduction can also be seen as a perturbation of the (2, 0)-phase reduction.

Dynamics
In this section we consider two different orbits, i.e., the synchronized orbit and the splay orbit, and compare their stability in a few different systems, including the full system (3.3) and various phase reductions on different order.

Synchronized State
First, we consider the synchronized state in the full system.A state in the phase space is called synchronized if all the oscillators are at the same position.In the full system (3.3) this state is defined by Consequently, if a state is synchronized, it is uniquely given by its amplitude R ⋆ := R i and its phase φ ⋆ := φ i for any i = 1, . . ., N .Due to the S N symmetry of (3.3) the set of synchronized states (4.1) is dynamically invariant.Thus, we can insert the ansatz (4.1) into the system (3.3) to obtain ODEs for R ⋆ and φ ⋆ as Based on these equations one can see that always R ⋆ = 1 on the invariant torus and that the rate of attraction to R ⋆ = 1 is given by m(1 + δg(φ ⋆ )) 2 .Since |δ| is small, this rate is mostly governed by m < 0.Moreover, the phase φ ⋆ evolves with constant speed ω.Consequently, the synchronized orbit on the invariant torus is given by γ f (t) = (½, (ω 0 + ωt)½), which is periodic with period T = 2π/ω.Now, we look at synchronized states in phase-reduced systems.In a phase-reduced system, there are no amplitudes and thus, a synchronized state is present when the single condition Here, a synchronized state is only determined by its phase φ ⋆ := φ i for any i = 1, . . ., N .Similarly to the full system, phase-reduced systems retain the S N symmetry and therefore the set of synchronized states is dynamically invariant.Inserting φ ≡ φ ⋆ into any of the phase reductions derived in Section 3 yields φ⋆ = ω.
Having established representations for the synchronized orbits, we now investigate their stability.Usually, when one wants to check the stability of a synchronized orbit, one changes into a co-rotating system, in which each synchronized state is an equilibrium.Then, one linearizes the vector field around this equilibrium and calculates the eigenvalues of this linearization.If they are all negative, apart from a single 0 eigenvalue that corresponds to perturbations along the continuum of synchronized states, the synchronized orbit is linearly stable and linearly (neutrally) unstable otherwise.However, when δ = 0, we cannot change into a co-rotating coordinate system, since the full system as well as phase-reduced systems are not S symmetric.In particular, the rate of attraction to the limit cycle depends on the position on the limit cycle.Consequently, one needs to take averages over all rates of attraction of one period of the limit cycle.These averages are called Floquet exponents.The concept of Floquet exponents and Poincaré return maps is often helpful when analyzing the stability of periodic orbits [15,39].To understand it, let us consider the general differential equation where H is a smooth vector field on the phase space X and TX is the tangent bundle.Suppose there is a periodic orbit γ : [0, T ] → X with period T .Later we want to apply this to the full system, in which x = (R, φ) T ∈ R N ≥0 × S N , and phase-reduced systems with x = φ ∈ S N .To analyze the stability of the orbit γ, one assumes a perturbation of the starting point of the periodic orbit x(0) = γ(0) + ǫη(0) and continues this perturbation along the periodic orbit such that x(t) = γ(t) + ǫη(t).One the one hand, if all possible perturbations η(0) decay after one period T , we expect the periodic orbit γ to be stable.On the other hand, if some perturbations η(0) grow in amplitude, the orbit γ is unstable.Therefore, we solve for η(t).However, since that is difficult to do in general, we first linearize in ǫ to obtain on first order which is a system of linear ODEs with a time-dependent coefficient matrix A(t).Now, let Φ(t) ∈ R N ×N be a fundamental solution of this ODE such that η(t) = Φ(t)η(0).To determine the linear stability of the periodic orbit γ(t) one propagates all possible perturbations η(0) over one period T of the orbit and then looks at the eigenvalues of the map Φ(T ).Of course, this map has an eigenvalue 1, that corresponds to the eigenvector that represents a perturbation along the periodic orbit.All other eigenvalues λ k , k = 1, . . ., N − 1 are the multipliers of a Poincaré return map.They are related to the Floquet exponents q k of the orbit as λ k = e T q k , for k = 1, . . ., N − 1.If the largest absolute value of all Poincaré return map multipliers (PRMMs) is less than 1, the periodic orbit is linearly stable.If one PRMM has an absolute value greater than 1, the orbit is linearly unstable.Next, we apply this concept to the synchronized orbit in phase-reduced systems.We start by calculating the PRMMs for the (1, ∞)-phase-reduced system (3.6).When putting this system into the framework of (4.2), we see that the matrix A(t) is where ½ is a N × N matrix where all entries are ones and I is the N × N dimensional identity matrix.Due to the S N symmetry of the synchronized state in the phase-reduced systems, the matrix A(t) has the special property that it is just a multiple of ½ − N I. Since matrices of this form commute with each other, a fundamental solution Φ(t) of (4.2) can explicitly be calculated using the matrix exponential: Integrating f (γ(t)) over one period of the orbit γ(t) = ω 0 + ωt yields Combining this with the fact that the eigenvalues of exp( a N (½ − N I)) are e −a with multiplicity N − 1 and 1 with multiplicity 1, we infer that the critical PRMM is Interestingly, this is independent of δ even though the system (1, ∞)-phase reduction includes all orders in δ.Consequently, the stability of a synchronized orbit is unaffected by δ in a phase reduction without higher-order interactions.A similar calculation yields that the critical PRMM of the synchronized orbit in a (2, 0)-phase reduction is The critical PRMM in a (2, 1)-phase reduction agrees with (4.3), which can be shown using the formulas (3.13) and (3.15b) if g(φ) = sin(φ).The independence of λ crit on δ in a (2, 1)phase reduction can be explained by a symmetry mapping δ to −δ.Assuming g(φ) = sin(φ), the symmetry would alter the system such that the limit cycle is then parameterized by r = 1−δg(φ).
Using the S-symmetry of the system one can apply a phase shift φ → φ + π to obtain the original system where the limit cycle is given by r = 1 − δg(φ + π) = 1 + δg(φ).Consequently, the S-symmetry causes a Z 2 symmetry mapping δ to −δ.Contributions of δ to λ crit can thus only be via even powers of δ, which is why it δ does not appear in the critical PRMM of a (2, 1)-phase reduction.When g is a general function, λ crit of a (2, 1)-phase reduction still does not depend on δ and the formulas to show that are found in the Appendix.The stability of the synchronized orbit in a (2, 1)-phase reduction thus agrees with its stability in a (2, 0)-phase reduction, i.e., one for which the limit cycle is circular.Thus, when concerned with the stability of the synchronized orbit, deformations of the limit cycle can be ignored to first order.Finally, when g(φ) = sin(φ), the critical PRMM in a (2, 2)-phase reduction it is given by which finally shows the dependence on δ.
Having derived stability conditions for synchronized orbits in phase-reduced system, we now analyze the stability of the synchronized orbit in the full system.Unfortunately, when applying the concept (4.2) to the full system, the matrices A(t) and A(s) do not commute with each other.Therefore, it is not possible to use the matrix exponential to analytically compute PRMMs or Floquet exponents, but we need to use numerical methods to determine them.Yet, in the special case δ = 0, the stability analysis of these periodic orbits simplifies quite significantly.In fact, in this case, the full system has an S symmetry, which acts by shifting all oscillators by a constant phase.Then, one can also change to a co-rotating coordinate frame, in which ω = 0.In these new coordinates all synchronized states are then equilibria.The spectrum of the right-hand side of the system at the synchronized state then contains information about the stability.There will be one 0 eigenvalue, since there is a one-dimensional continuum of synchronized states.If all the other eigenvalues have negative real part, the synchronized state as an equilibrium in the co-rotating frame is linearly stable and thus the synchronized orbit in the original system inherits this stability.Conversely, if one eigenvalue has positive real part, the synchronized orbit in the original system is unstable.Conducting this analysis for the full system (3.3)yields that the linearization of the right-hand side is given by a matrix The eigenvalues of this matrix can explicitly be calculated and are given by where we denote them by q k for k = 1, . . ., 2N because they describe the instantaneous rate of attraction to the periodic orbit and thus relate to the Floquet exponents.In fact, since this instantaneous rate of attraction is constant over the whole orbit, these eigenvalues agree with the Floquet exponents.While the first N eigenvalues correspond to perturbations of the phases φ, the last N eigenvalues originate from perturbations in the radial directions.To compare this model with phase-reduced models, we assume that −m is big enough such that the last N eigenvalues can be neglected and the critical Floquet exponent q crit is given by q 2 , . . ., q N .Of course there is also the zero eigenvalue q 1 .However, that does not contribute to the stability as it corresponds to a perturbation along the continuum of synchronized states.Given the critical Floquet exponent, one can then obtain the critical PRMM by simply calculating λ crit = exp(q crit T ) = exp( 2π ω q crit ).When δ = 0 in the full system, PRMMs can only numerically be calculated, as shown in Figure 2(a).Figures 2(b-d) compare PRMMs from the full system with PRMMs from phase-reduced systems.

Splay State
While all oscillators gather at one point on the circle if they are synchronized, one can say that the splay state is the opposite of that.A splay state is given when the oscillators phases are equidistantly distributed on the circle.More specifically, a state φ ∈ S N is a splay state if there is a permutation σ : for k = 1, . . ., N − 1.By relabeling the nodes, we might also assume that σ is the identity map.Thus the set of splay states is given by

.5)
A splay state in this set can be characterized by just the first phase φ 1 .This set has a Z N := Z/(N Z) symmetry group that acts on a state by shifting the indices, which have to be understood modulo N , of each oscillator by a constant integer.Now, let us consider the set of splay states in phase-reduced systems with δ = 0. Since the right-hand sides of (1, 0)-and (2, 0)-phasereduced systems are equivariant with respect to this group action, the set of splay states is dynamically invariant.In particular, when inserting a splay state into the right-hand side of a (1, 0)-phase reduction and a (2, 0)-phase reduction, it follows that φk = ω − K sin(α) =: ω.Therefore, if ω = 0, there exists a periodic orbit γ(t) ∈ S N with γ k (t) := ωt + 2πk N , that has period T = 2π ω = 2π ω−K sin(α) .We refer to this orbit as the splay orbit.Next, we consider splay states in the full system (3.3) and still assume δ = 0. Inserting the ansatz (4.5) into the full system (3.3)yields that the amplitudes on the invariant torus are given by Therefore, we call a state (R, φ) ∈ R N ≥0 ×S N a splay state if φ ∈ D and the amplitudes satisfy (4.6).Then, the splay state in the full system has the same symmetry group Z N as splay states in phasereduced systems.Moreover, since the right-hand side of the full system (3.3) is again equivariant with respect to this symmetry group, the splay state is dynamically invariant.Furthermore, the angular frequency of the phases is φk = ω−K sin(α) = ω as for phase-reduced systems.Therefore, there there exists a periodic orbit γ(t) = (R(t), φ(t)) T with R k (t) = R ⋆ and φ k (t) = ωt + 2πk N , that has the same period as the one in phase-reduced systems.
When analyzing the stability of these splay orbits in both phase-reduced and the full system, it is important to note that splay orbits are just one single periodic orbit in a whole continuum of periodic orbits.In particular, in phase-reduced systems, all incoherent states, that are characterized by rotate around the circle with constant frequency ω and are a union of manifolds [2].In the full system states (R, φ) ∈ R N ≥0 × S N with R k = R ⋆ for k = 1, . . ., N and Z = 0 rotate around the circle with the same constant frequency.Therefore, there exist further periodic orbits, which we refer to as incoherent orbits.Since a splay state is a special incoherent state, but in general the set of incoherent states is larger than the set of splay states, the splay orbit is only one orbit in a continuum of incoherent orbits.Consequently, when analyzing the stability of the splay orbits with PRMMs or Floquet exponents, there is always one neutral multiplier or exponent, respectively.The only exception is present when the set of incoherent states coincides with the set of splay states, i.e., when N = 3.To overcome the problem of neutral stability, we restrict ourselves to N = 3.
If δ = 0, we can change into a rotating frame coordinate system, in which splay states are equilibria.After having done that, we linearize the right-hand side at the splay state and thereby obtain Jacobians with eigenvalues in a (1, 0)-phase reduction and Ke ±iα in a (2, 0)-phase reduction.An analytical derivation of the eigenvalues in the full system turned out to be too complicated.In both cases, the critical Floquet exponent is given by q crit = q 2,3 and thus the critical PRMM is λ crit = e T q crit .While all the previous theory was only valid for δ = 0, we now investigate what happens if δ = 0.In this case, the right-hand sides of both the phase-reduced systems and the full system are no longer equivariant with respect to Z N .Therefore, splay states are in general not invariant anymore.However, when N = 3 and δ = 0, there is a single periodic orbit, i.e., the splay orbit.For general parameter values, this orbit has no Floquet exponents whose absolute value equals 1.Thus, this splay orbit is hyperbolic.Slight changes in δ away from 0 preserve the existence of a periodic orbit in the neighborhood of the splay orbit.To illustrate the stability of these orbits, we numerically search for periodic orbits in a neighborhood of the splay orbit in the full system and phase-reduced systems.Then, we numerically calculate their critical PRMMs to determine the stability of these orbits, see Figure 3.
A numerical analysis revealed that there is a subcritical Neimark-Sacker bifurcation in Figure 3(f) when K is negative and the modulus of the critical PRMM passes through 1.In particular, there are two complex conjugated PRMMs that pass through the complex unit circle.This correctly represents the bifurcation behavior of the full system, as in Figure 3(d).A (1, ∞)-phase reduction does not even capture the bifurcation, see Figure 3(e).The bifurcation at K = 0 is degenerate.When K = 0 there is no coupling and so all eigenvalues are 0.

Phase Reduction Beyond All-To-All Coupled Networks
In Sections 3 and 4, we have started with a system of coupled oscillators, derived various phase reductions and compared the stability of synchronized and splay orbits.The system we started with (3.2), consists of N complex Stuart-Landau oscillators coupled with each other via a meanfield coupling, i.e., each oscillator influences every other oscillator in the same way.However, instead of this all-to-all coupling, one could also assume that the coupling between oscillators is described by a graph.In Section 5.1 we adopt our calculations from the previous sections to a non all-to-all coupling and Section 5.2 then contains an interpretation of the higher-order interactions terms that we derive from a phase reduction of nonlocally coupled Stuart-Landau oscillators.This provides additional evidence that it is important to study dynamical systems on hypergraphs which have recently been encountered in many different contexts; see for example [38,12,37,23,24,20].

Derivation of the phase-reduced Dynamics
Now, let us assume that the coupling structure is not given by an all-to-all topology, but instead there by a (possibly directed and weighted) graph Γ = (V, E) that is described by its adjacency matrix A ∈ R N ×N with entries a kl .Then, the governing equation is which contains (3.2) as a special case when a kl = 1 for all k, l = 1, . . ., N .Now, one can do the same analysis as in Section 3, but with (5.1) as a starting point.The procedure from Section 3 is directly applicable to (5.1),only the resulting formulas are slightly different.Therefore, will not explain the whole procedure again, but only state the results.
The system (5.1) can be written in polar coordinates A k = r k e iφ k .Transforming the radii as and that gradients on higher order in δ are generally of the form . . .

Second-Order Phase Reductions as Higher-Order Networks
We now discuss the individual coupling terms that constitute the (2, 0)-phase reduction.The coupling includes nonpairwise terms and we discuss how the coupling can be interpreted as a higher-order phase oscillator network on hypergraphs that can be derived from the original graph Γ = (V, E) that describes the coupling of the nonlinear oscillators.In summary, the (2, 0)-phase reduction is given by where agrees with (5.3) for δ = 0 and P (2,0) k (φ) as specified in (5.4) contains the second order terms.First, note that the coupling of the first-order phase reduction (5.5) is posed on the graph Γ (1) := Γ that describes the interactions of the coupled nonlinear oscillator network.
The second-order phase interaction terms (5.4) not only contain pairwise interactions along network edges but also nonpairwise interactions between triplets of oscillators.Collecting all terms, the phase interactions can be interpreted as interactions on two 3-uniform directed hypergraph with the (adjacency) 3-tensors ĥ, h ∈ R N ×N ×N with coefficients ĥkli := a kl a ki , hkli := a kl a li , where a triplet (k; l, i) corresponds to a directed hyperedge with tail k and head {l, i}.
3 The expression for P (φ) can be generated with the Mathematica code accompanying this paper, see [11].
Figure 4: Illustration of second order phase interactions in K that affect node k (indicated in red).The first row corresponds to terms appearing in ĝ, while the second row lists terms of ḡ.The black lines indicate edges in the graph Γ that describe the interaction of the unreduced nonlinear oscillator network.The (hyper)edges that describe the directed phase interactions are indicated by yellow blobs (with node k being the head).These include pairwise interactions (panels a1, a2), pairwise interactions that may be virtual (panel b2), and three types of nonpairwise interactions that describe the nonlinear influence of nodes i, l onto node k: One does not depend on the state of node k itself (panel b1) and two that depend on the phases of all three nodes (panels c1, c2).
The hypergraph capturing nonpairwise higher-order interactions 'inherits' directionality from the original graph Γ: The coefficient h kli is nonzero if l, i are both in the neighborhood of the target node k in Γ (cf. Figure 4(1)), while hkli is nonzero if there is a path from i via l to k in Γ (cf. Figure 4(2)).The coupling functions along these hyperedges are (5.6) These can be obtained from (5.4) using trigonometric identities and we have Since the coupling functions (5.6) contain both pairwise and nonpairwise phase interactions, the interaction structure can be broken down term separating pairwise and nonpairwise interactions: Each of the three terms can be associated with a particular type of interaction, which results in six subclasses in total shown in Figure 4.
First, there are pairwise correction terms to the first-order phase reduction, that correspond to the first term in the definition of ĝ and ḡ; see (5.6).The coupling of these pairwise correction terms is posed on the graph Γ (2) a := Γ that is the same graph as the coupling of the full system and the coupling of the first-order phase reduction Γ (1) ; see Figure 4(a).One of these pairwise correction terms is weighted with the degree of the node k; see Figure 4(a1), while the other term is weighted with the degree of node l (that is in the neighborhood of k), see Figure 4(a2).Moreover, the interaction function sin(φ l − φ k + α) agrees with the interaction function of a firstorder phase reduction (5.5) up to the shift sin(α) and a constant factor.Since the two pairwise correction terms in Figure 4(a) share the same coupling structure Γ (2) a and the same interaction function, they can be combined into Based on this formula, one can see that the sign of the pairwise correction terms is determined by comparing the degree of node k with the average degrees of all neighbors l of k.
Second, consider the second term in ḡ that describes a pairwise interaction between node k and node i that is in the neighborhood of a the neighbor l of k.While this yields a pairwise interaction from node i to node k, there may not necessarily be and edge (i, k) ∈ E(Γ) from i to k in the graph Γ that describes the original network of coupled nonlinear oscillators.If (i, k) ∈ E(Γ) then this second-order term describes a second-order correction to the first-order interaction.If (i, k) ∈ E(Γ) then this interaction can be considered as a virtual edge, which is present in a weighted graph Γ (2) b2 defined by the adjacency matrix C = (c ki ) k,i=1,...,N with coefficients c ki := l a kl a li but not in Γ.As one can see from the definition of the entries c ki of the adjacency matrix of Γ (2) b2 , this interaction is weighted by the number of paths connecting k to i in the graph Γ.
Third, the second term of ĝ represents a triplet interaction where nodes i and l-both neighbors of k-influence the node k; see Figure 4(b1).While the interaction is a nonpairwise interaction involving three distinct nodes (i, l and the target k), this may be considered a 'nonstandard' nonpairwise interaction as the coupling is independent of the state of node k.This interaction can be represented by a directed and possibly weighted hypergraph H (2) b1 represented by a corresponding adjacency 3-tensor indexed by k, l, i.Note that there will be a symmetry that allows swapping the indices l and i, but in general this hypergraph is still directed as one can not arbitrarily permute all indices k, l, i.
Finally, there are two further triplet interactions.The first triplet interaction, i.e., the last summand in the definition of ĝ is an interaction between two neighbors i and l of the node k.The coupling structure of this interaction can be described by a directed and weighted hypergraph H (2) c1 , whose 3-tensor agrees with ĥ.Note that there is a symmetry between i and l, but one cannot arbitrarily permute all indices which is why in general the hypergraph is directed.The second triplet interaction is governed by the last summand in the definition of ḡ, is one between the node k itself, a neighbor l of k and a neighbor i of l.Again, the coupling structure can be described by a directed hypergraph H (2) c2 .This time, however, the 3-tensor that describes the hypergraph corresponds with h and in general, there this hypergraph does not possess any symmetry with respect to a permutation of indices.
To summarize, the first-order phase reduction, that we have considered, consists of only one interaction term, whose coupling structure Γ (1) agrees with the coupling structure Γ of the full system (5.2).In a second order phase reduction, quite a few new interaction terms appear.While the coupling of some of them is posed on a graph that agrees with Γ, the coupling of others is determined by a graph that consists of virtual edges that might not be present in Γ.Moreover, there are also three types of triplet interactions on directed hypergraphs, which can be derived from the adjacency matrix of Γ.
To conclude this section, we want to remark that a second order phase reduction contains interaction terms on directed hypergraphs, even when the underlying graph Γ, that determined the coupling in the full system (5.2), is undirected and unweighted.The only exception is when Γ itself is an all-to-all graph.However, whenever Γ is connected, yet non all-to-all, there exists an open triangle as seen in Figure 4(c1), which causes the hypergraph, that governs the second order phase reduction, to be directed.

Discussion
Phase reductions provide a useful tool to analyze the dynamics of coupled oscillator networks.Here we derived explicit expressions for nonlinear oscillations with phase-dependent amplitude subject to simple diffusive coupling.By using a suitable coordinate transformation, our results also apply to systems where the limit cycle is simple but the coupling is phase dependent or a combination thereof.The reduced phase equations allow to analyze, for example, phase dynamics for oscillators that are-if uncoupled-further away from a Hopf bifurcation point where the limit cycle becomes noncircular in general.
While the shape of the limit cycle affects the collective dynamics, a first-order phase reduction is insufficient to capture the dynamical effects of the amplitude dependence: The phase reduction needs to be at least of second order in both the coupling strength K and the parameter δ that describes the perturbation from a circular limit cycle.We showed that second-order phase reductions were able to accurately predict the stability properties of the synchronized and splay orbit when all terms of up to second order in K and δ were included.Importantly, the amplitude dependence breaks the rotational symmetry of phase equations that is typical, for example, of the Kuramoto equations.While such symmetry breaking has been analyzed from the perspective of the phase equations [13], in our setup it arises through a perturbation of the underlying nonlinear oscillator.This perspective allows to make direct comparisons between the nonlinear system and the phase-reduced dynamics in contrast to phase reductions via normal forms [5], where normal form symmetries-that may be absent in the full equations [16]-can appear in the phase reduction.
As detailed in Section 5.2, the second-order phase interaction term for coupled oscillators on a given graph Γ can be interpreted as phase oscillator dynamics on (directed) hypergraphs.Specifically, the second-order phase interaction terms correspond to second-order corrections to interactions along edges of Γ, possible virtual pairwise connections between oscillator pairs that are not joined by an edge in Γ, and nonpairwise triplet interactions of different type.This highlights the importance to consider the dynamics on directed hypergraphs-these may appear implicitly (e.g., [7]) while explicit frameworks with distinct perspectives have been introduced in [1,19,41].Indeed, directionality of higher-order interactions has implications for synchrony [19] or the emergence of more complicated dynamical phenomena such as heteroclinic dynamics [7,10].
In principle, the analysis presented in Section 3 can be extended to derive higher-order interactions beyond second order in both the coupling strength K and the deviation parameter δ.Such higher-order interactions would include non-additive interactions between quadruplets of phases and allow to describe the approximate phase dynamics beyond second order.The main obstacle that has to be overcome is the algebraic complexity of these phase reductions.Already at second order, the terms for triplet interactions become quite complex that bring symbolic computer algebra software, such as Mathematica, to their limits.
Generalizing the computations to oscillator networks with nonlinear coupling, nonidentical oscillators, or strongly nonsinusoidal oscillations-all properties of real-world oscillator networkshighlights some of the challenges to compute phase reductions explicitly.First, nonlinear cou-pling makes equation (3.5) more challenging to solve (if an explicit solution is possible at all).If the nonlinear coupling contains nonlinear coupling between three oscillators (e.g., a coupling terms of the form A l A j A k )-corresponding to then we expect nonpairwise coupling already in the first-order phase reduction (see [4]).If the coupling is nonlinear in one oscillator, i.e., the coupling only includes coupling terms of the type A m j (as, for example, in [28]) then we have higher harmonics in (3.5) making it less tractable.Second, real-world oscillators are rarely perfectly identical motivating the question how heterogeneity affects the phase reduction.Here we assumed that the oscillators are identical and, in particular, that ω in (3.3b) is independent of k.
If ω depended on k the first-order phase reduction would not change at all and one could just replace ω by ω k everywhere.A second-order phase reduction, could theoretically also be derived by adapting the methods from Section 3. Practically, however, the terms of second order start to depend nonlinearly on ω = (ω 1 , . . ., ω N ).Thus, already when δ = 0 and N is small and finding a general solution for R (1,0) k (φ) is challenging.Moreover, due to the growing complexity of the second-order phase reduction, it is practically intractable.Extending this to δ > 0 only worsens the problem.A possible approach to overcome this problem is to assume that the intrinsic frequencies ω k are sampled from a probability distribution and consider a mean-field limit.Third, strongly nonlinear oscillations, such as FitzHugh-Nagumo neurons, are far away from a circular limit cycle (and would this require a large perturbation parameter δ).While this suggests the importance of expanding to high-order in δ, more direct approaches may be more appropriate [26] that allow to get a more qualitative understanding of the dynamics [3].

Figure 2 :
Figure 2: Comparison of the critical multipliers of the Poincaré return map of the synchronized orbit in the full system and phase-reduced systems.(a) shows the critical PRMM of the synchronized orbit in the full system in dependence of δ and the coupling strength K. Part (b) depicts the critical PRMM of the synchronized orbit in a (1, ∞)-phase reduction and part (c) illustrates the critical PRMM of the same orbit in a (2, 2)-phase reduction.Finally (d) depicts the critical PRMM in the full system, the (1, ∞)-phase reduction and the (2, 2)-phase reduction when δ = 0. Parameter values: ω = 1, m = −1, α = π/2 + 1/20, g(φ) = sin(φ).

Figure 3 :
Figure 3: Numerical calculation of periods and critical multipliers of the Poincaré return map of periodic orbits in a neighborhood of the splay orbit.The first column represents the full system, the second column is the (1, ∞)-phase reduction and the third column displays the (2, 2)-phase reduction.The upper row is the period of the resulting periodic orbit, while the lower row depicts the critical PRMM of this orbit.Parameters: α = π/2 + 1/20, m = −1, ω = 1, N = 3, g(φ) = sin(φ) ) and P