Identical phase oscillator networks: bifurcations, symmetry and reversibility for generalized coupling

For a system of coupled identical phase oscillators with full permutation symmetry, any broken symmetries in dynamical behaviour must come from spontaneous symmetry breaking, i.e. from the nonlinear dynamics of the system. The dynamics of phase differences for such a system depends only on the coupling (phase interaction) function $g(\varphi)$ and the number of oscillators $N$. This paper briefly reviews some results for such systems in the case of general coupling $g$ before exploring two cases in detail: (a) general two harmonic form: $g(\varphi)=q\sin(\varphi-\alpha)+r\sin(2\varphi-\beta)$ and $N$ small (b) the coupling $g$ is odd or even. We extend previously published bifurcation analyses to the general two harmonic case, and show for even $g$ that the dynamics of phase differences has a number of time-reversal symmetries. For the case of even $g$ with one harmonic it is known the system has $N-2$ constants of the motion. This is true for $N=4$ and any $g$, while for $N=4$ and more than two harmonics in $g$, we show the system must have fewer independent constants of the motion.


INTRODUCTION
The last 30 years has seen a lot of progress in understanding the collective behavior of coupled oscillating dynamical systems, especially since the pioneering work of Winfree [1] and the coupled phase oscillators of Kuramoto [2,3] which has inspired many studies of criteria for synchronization and applications in a variety of areas, notably neuroscience [4]; see also Acebrón et al. [5] and Strogatz [6] for reviews. Most of the theoretical work has focussed on the continuum limit of large populations of various types of oscillator, and has been guided by large numerical simulations-see for example Ullner and Politi [7].
Nonetheless, relatively small networks may show a lot of nontrivial and remarkable behavior that will be present in larger networks. We consider a general system of N coupled identical phase oscillators as in Kuramoto [2] but with a general coupling/phase interaction function g(ϕ): where (θ 1 , . . . , θ N ) ∈ T N represents the phase of the oscillators, ω the natural frequency of the oscillators and g(ϕ) is the 2π-periodic smooth coupling function (also called the phase interaction function). Such systems arise naturally as approximations of coupled limit cycle oscillators in a weak coupling limit; see for example Ashwin and Rodrigues [8].

Identical Coupled Oscillators and Symmetric Dynamics
This case of identical oscillators is clearly highly symmetricindeed, it is equivariant under the group of permutation symmetries S N acting by permutation of coordinates as well as being equivariant under translation along the diagonal (1, . . . , 1)-giving a symmetry group Ŵ = S N × T. Some consequences of these symmetries have been studied in Ashwin and Swift [9] and subsequently in Brown et al. [10]-tools from symmetric dynamics [11] allow one to understand a lot of the behavior of the system by looking at the action of the symmetry group on phase space. In particular, for any point in θ ∈ T N we associate the isotropy subgroup which is the subgroup g ∈ Ŵ gθ = θ ⊂ Ŵ. For each subgroup G ⊂ Ŵ we can define the fixed point subspace Fix(G) = θ gθ = θ ∀ g ∈ G ; it is a standard result from Golubitsky et al. [11] that for any G, Fix(G) is dynamically invariant for any Ŵ-equivariant dynamical system. There is a twist in the story in that here Ŵ acts on the manifold T N rather than R N , giving rise to more non-trivial structures than expected from Golubitsky et al. [11]. In this paper we also explore the presence of time-reversal symmetries [12,13] for even g.
In order to understand how (Equation 1) may have complex spatio-temporal attractors of various symmetry types, we study this fully symmetric case for small N and specific g to gain an insight into the general case. For the case of Kuramoto-Sakaguchi coupling [14] g(ϕ) = − sin(ϕ − α) (2) one can perform a detailed analysis of the general N case and one finds either synchrony or antisynchrony are attracting, with the cases α = π/2, 3π/2 giving special integrable structure [15] that corresponds to "vertical branches" at bifurcation. However, this sort of behavior is highly degenerate for a general dissipative system. On addition of higher harmonics to the coupling function we expect, and usually find, that the degeneracies typically unfold into generic situations. Indeed such higher harmonic terms are present in a generic system even near Hopf bifurcation [8,16]. A range of solutions, their stabilities and bifurcations, and the invariant subspaces imposed on phase space for the case of general g is discussed in Ashwin and Swift [9]. This paper briefly reviews and then extends and applies some of these results in two contexts. The first aim of the paper is to consider the detailed bifurcation behavior for the most generic two harmonic coupling function. This means we determine which degeneracies of Equation (2) are unfolded by addition of the second harmonic, for the cases of N ≤ 4. In doing so we extend the analysis in Ashwin et al. [17] to the most general two harmonic coupling. The paper [17] examines the cases N = 3 and 4 for Equation (1) and presents complete global bifurcation diagrams for the phase interaction function discussed in Hansel et al. [18] g(ϕ) = − sin(ϕ + α) + r sin 2ϕ. (3) However, for this parameterized family the only even coupling functions are already in Equation (2). To next order [8], the generic lowest order terms will be of the form where q, r, α and β are arbitrary constants (specific cases of such coupling was considered by [19] in the context of N = 5). More generally one can consider as in Daido [20] L-harmonic coupling where, without loss of generality, we can assume q 1 q L = 0, and Equation (4) corresponds to L = 2.
The second aim of this paper is to consider the behavior of the system (Equation 1) for arbitrary, but even g(ϕ) = g(−ϕ). Following Bick [21], we note there will be special structure in the dynamics, and a range of non-trivial recurrent dynamics appears for g that are small perturbations of even functions. This has been illustrated by the observation of chaotic attractors [22] and extreme sensitivity to detuning [17] for fully symmetric systems of oscillators. Indeed, chimera states [23][24][25] for structured (but not necessarily fully symmetric) networks are usually found close to even coupling. We observe that there is zero divergence for Equation (1) and time-reversal symmetries for the system of phase difference (though typically not for the original system). It is known there are N − 2 independent smooth integrals of the motion for any N ≥ 3 for the special case g(ϕ) = cos(ϕ) [15]. This still holds for general even g and N = 3 but, surprisingly, we show here that it does not hold for N = 4 and general even g: the flow is not fully integrable.
The paper is structured as follows: after discussing some parameter and phase symmetries, we recall the reduction of the dynamics to a canonical invariant region. In Section 2 we summarize some stability properties and bifurcations of periodic orbits by symmetry type. Section 3 applies this to general two harmonic coupling and N = 3 or 4 oscillators, extending the bifurcation analysis in Ashwin et al. [17] to general second harmonic coupling. This includes a new characterization of bifurcation of two-cluster states as the discriminant of a certain polynomial equation. Section 4 examines some consequences of the coupling function g itself being odd, or even for small N. If g is odd, it is known that the system for phase differences (but not the original system) has a variational structure and can be written as a gradient flow. If g is even, there are additional timereversal symmetries [13] for the system for phase differences, as well as a zero divergence condition. We characterize the points with nontrivial time-reversal symmetries for N = 3 and N = 4 in detail. Finally, Section 5 considers the relation between even coupling and integrability. For N = 3 and general even g, or for N ≥ 4 and g(ϕ) = cos(ϕ) it is known since [15] that there are N − 2 independent integrals of the motion. This constrains the dynamics for N ≥ 4 to be at most two-frequency quasiperiodic. We show this is not the case for N ≥ 4 and more general even g: for examples with N = 4 and g with L ≤ 4 harmonics, we find structures that respect the time-reversal symmetries, but that cannot appear in a fully integrable structure. We discuss this, and highlight some open questions.
Using Equation (6) and rescaling time by a positive scaling, we assume from hereon that q = −1: using Equation (7) we can assume from hereon that r ≥ 0. We choose q = −1 to give easy comparability to Ashwin et al. [17]. There is a timereversal parameter symmetry (such that g(ϕ) → −g(ϕ) for all ϕ), assuming q = −1 this is (r, α, β) = (r, α + π, β + π) By classifying qualitative dynamics up to time-reversal, we will can reduce to parameters in the region in the later bifurcation analyses; using Equation (8) we can classify for β ∈ [π, 2π). Finally, there is a phase-reversal parameter symmetry (i.e., such that g(ϕ) → −g(−ϕ) for all ϕ), assuming q = −1 this is This means in principle we can reduce to Equation (9) with β ∈ [0, π/2] and apply phase and time reversal as necessary.

Phase Symmetries: the Canonical Invariant Region
Because the dynamics of Equation (1) depends only on phase differences, one can usefully reduce the dynamics from a flow on T N to one on T N−1 which represents only the phase differences. This corresponds to looking at the quotient dynamics of the system on the T-orbits. More precisely, if we consider a projection : T N → T N−1 that maps the T-orbits onto points and considerθ = (θ ) ∈ T N−1 , the original system (Equation 1) becomes a reduced system on T N−1 of the form This clearly depends on g and N but not on ω.
There are several ways to choose , or, equivalently, a set of N − 1 phase differences θ j − θ k that span the set of all phase differences. However, permutation symmetry will never act orthogonally on such a set of phase differences if N ≥ 3, meaning that symmetries are not obvious from plots of phase differences. For N = 3 and N = 4 there is a useful representation introduced by Ashwin and Swift [9] that avoids this problem. If we examine an orthogonal projection of the dynamics onto a codimension one subspace orthogonal to the diagonal (1, . . . , 1), this will carry an orthogonal (and irreducible) representation of S N on R N−1 , the covering space of T N−1 . In addition, the complement of fixed point subspaces where two of the phases are identical (isotropy S 2 ) forms a partition of R N−1 into connected components that are all images of the canonical invariant region (or CIR) [9]. This is the set whose boundary consists of points with S 2 isotropy. There is a residual symmetry Z N = Z/NZ = τ on C generated by Note that splay = (0, θ 2 , . . . , θ N ) ∈ C with θ k = θ k − 1 + 2π N is the only fixed point of the action of τ in C [9]. Figure 1 illustrates the C for the cases N = 3 and N = 4the boundaries of C are invariant for dynamics governed by Equation (1) for any coupling function g, though they may not be invariant for arbitrarily small perturbations that make the g or the ω different for different oscillators-in this case we may have extreme sensitivity to detuning [17]. Nonetheless, even in this case the trajectory will remain most of the time in one image of the CIR, and the dynamics can usefully be understood as being a perturbation of the highly symmetric case.
Within C there are representatives of points with all possible symmetry types. Particularly important are ℓ-cluster states which correspond to states where for some partition of N = m 1 + · · · + m ℓ there are precisely m k oscillators at the same phase, for some ℓ ≥ 2 (the number of clusters) and m k ≥ 1 (the size of the kth cluster). These correspond to fixed point subspaces for isotropy groups S m 1 × · · · × S m ℓ . These are located on the boundary of C-there are also points with spatio-temporal symmetries (S m 1 × · · · × S m ℓ ) M × s Z M with M > 1 and N = M(m 1 + · · · + m ℓ ); see Ashwin and Swift [9]. NB: To simplify notation we omit copies of S 1 when describing cluster states where this is clear from context: for example S 2 for the case N = 4 corresponds to a three cluster state with isotropy group S 2 × S 1 × S 1 .

STABILITIES AND BIFURCATIONS OF SOLUTIONS FORCED BY SYMMETRIES
We recall from Ashwin and Swift [9] a number of important synchronous solutions that are forced by symmetry, and their bifurcation properties. In (C,D) the solid lines have isotropy S 3 × S 1 while the long dashes have isotropy S 2 × S 2 . The short dashed lines have isotropy (S 1 ) 2 × s Z 2 -typical points being (a, b, a + π, b + π ). In each case there is a residual Z N−1 symmetry indicated by the arrows in (B,D) and (N − 1)! symmetric copies of C pack a generating region for the torus.

Stability and Bifurcation of Synchronized Solutions
The fully synchronous solution (in-phase oscillation) of Equation (1) is a periodic orbit given by where θ (t) = sync t + θ 0 and sync = ω + g(0).

Lemma 1 ([9]). The fully synchronous solution sync (t) is linearly stable if and only if
For general two harmonic coupling (Equation 4) with q = −1, we have g(0) = 2r cos β − cos α and so there is bifurcation to loss of stability of this solution when independent of N.
A splay phase solution (rotating wave) is, up to permutation, a periodic orbit of the form where θ (t) = splay t + θ 0 . There are (N − 1)! different splay phase solutions from the possible permutations, but only one in the canonical invariant region C. Let i = √ −1.

Lemma 2 ([9]
). The splay phase solutions are linearly stable if the real parts of Let η j = g ′ 2π N j and ν p = ρ p = exp( 2π i N p). If = splay then the matrix given by Note that ν j N = exp 2π i N jN = 1 implies λ N ( ) = 0 for any splay phases. The system has [N/2] complex conjugate pairs The system for phase differences near the splay state has an Andronov-Hopf bifurcation (which can be degenerate or multiple) when Re(λ p ( )) = 0 for some index p = 1, . . . N − 1, p = N/2, and Im(λ p ( )) = 0. For general two harmonic coupling (Equation 4) this can be expressed as conditions on the parameters q, r, α, β, though note that unlike the synchronous state, this will depend on N in a fairly complex way. The cases for N = 3 and N = 4 are discussed in Section 3.

Antiphase States: the Set with Zero Order Parameter.
We discuss briefly a set that is a union of invariant manifolds for some forms of coupling or for small enough N. We note that it is not invariant for generalized coupling N ≥ 5. Consider the Kuramoto order parameter for system (Equation 1) is and define the set with zero order parameter R = 0 by ; it can be characterized by the set where two variables are given in terms of the other N − 2 variables by expressions The manifolds can be written in phase differences (for example However, although some solutions (in particular, splay states Equation 17) belong to M (N) for arbitrary coupling function and arbitrary g, quite surprisingly it is not generally invariant for general coupling. Proof. For N = 2, the equation exp(iθ 1 ) + exp(iθ 2 ) = 0 has solution θ 2 = θ 1 + π, while for N = 3, the equation Both cases correspond to the splay phase state which is a zero dimensional fixed point subspace for the phase differenceshence it is dynamically invariant. For N = 4, the equation This corresponds to the fixed point subspace of (S 1 × S 1 ) 2 × s Z 2 which is also therefore invariant.
On the other hand, in the case N ≥ 5 the corresponding sets M (N) are not fixed point subspaces, and therefore not necessarily invariant. For example, the set M (5) includes points of the form which is a two dimensional subspace on which However, the smallest fixed point subspaces the contains all points of this form is the trivial subspace T 5 , meaning the symmetry does not force invariance of M (5) .
Note that for N = 4 coupled oscillators, the set (one half-line in each invariant region) as a union of invariant manifolds that correspond to the states with isotropy Z 2 .

Stability and Bifurcation of Two Cluster States
Consider any two-cluster states with isotropy S p × S N − p given by such that θ i = θ a (t) for i = 2, . . . , p and θ i = θ b (t) for i = p + 1, . . . , N. If we denote ϕ p = θ a − θ b then (see for example Orosz et al. [27] for the general case) Equation (1) gives In this phase difference coordinate we can write Any fixed point in the phase differences will therefore satisfy For general two harmonic coupling (Equation 4) this means that Frontiers in Applied Mathematics and Statistics | www.frontiersin.org This has a trivial solution: the synchronous state with ϕ p = 0. Removing the factor of sin(ϕ p /2) this means the nontrivial two-cluster states will have a phase difference ϕ p where 0 = (N − 2p) sin(ϕ p /2) q sin α + 4r sin β cos 2 (ϕ p /2) + N cos(ϕ p /2) q cos α + 2r cos β 2 cos 2 (ϕ p /2) − 1 . (23) Note that two cluster-states for g(ϕ) = − sin(ϕ − α) (i.e., the special case r = 0) are considered in Burylko and Pikovsky [26][ Equation 17]. By setting τ = tan ϕ p /2, then Equation (23) can be expressed as a cubic equation. If we write p = NP so that 0 < P < 1, this implies the bifurcations occur precisely when the cubic equation for τ has discriminant zero, namely when where a k are polynomials in q, r, and cos and sin of α and β. The expressions for the a k are given in Appendix A (Supplementary Material).
There are also bifurcations with fixed points of twocluster manifold that occur in transversal to this manifold direction: these create three-cluster states. For example, pitchfork bifurcations (bifurcation line ID in Ashwin et al. [17], Figure 1 for three oscillators can occur transverse to the lines ϕ 1 = 0 (or ϕ 2 = 0, ϕ 1 = ϕ 2 ). Such bifurcations cannot be found by examination just of the equation on invariant manifold (Equation 20) or (Equation 21) but need to consider the transverse eigenvalues of the solutions, expressions for which are given in Orosz et al. [27]. Later on we use numerical path following to locate these.

BIFURCATION SCENARIOS FOR LOW N AND TWO HARMONIC COUPLING
We review and extend the exploration of parameter space of Equation (1) with general two harmonic coupling Equation (4) from the special case β = 0 [17] for N = 3 and N = 4 to more general β = 0. Figure 2 illustrates the bifurcations that determine the structure of the phase space for N = 3 in the case q = −1, on varying α ∈ [0, π) and r ≥ 0, for four choices of β. Note that the case β = 0 in the dot-dashed region corresponds to the case considered in Ashwin et al. [17]. Some of the curves are calculated analytically while the remainder of the curves on Figure 2 are obtained by numerical path following of bifurcations in the reduced system (Equation 11) using XPPAUT [28].

Bifurcations for N = 3
Note that from the previous section, we know that the inphase solution undergoes bifurcation at (Equation 16). The splay phase solution with N = 3 the stability is governed by the eigenvalues We summarize the bifurcations for β = 0 in Figure 2. One of the more surprising observations of the analysis above can be summarized as follows.

Corollary 1.
In order for the Hopf bifurcation of splay phase and the bifurcation of in-phase to occur at different points in parameter space for Equation (1) with N = 3, it is necessary to have a coupling function g with at least three harmonics. Note that the case β = 0 in the dot-dashed region corresponds to the case considered in Ashwin et al. [17] for N = 4. As before, the in-phase solution has a bifurcation at locations where (Equation 16) is satisfied. For the splay phase solution with N = 4 the stability transverse to the rotating block (0, π, θ, θ + π) is governed by the complex eigenvalues  Figure 3. The dot-dashed regions for the case β = 0 corresponds to the diagram detailed in Ashwin et al. [17]. Note that in the case β = π/2 and α = π/2 or 3π/2 the system is integrable, while the case β = 3π/4 can be found from that for π/4 by the phase-and time-reversal symmetry (r, α, β) → (r, π − α, π − β).

Bifurcations for N = 4
meaning there is a Hopf bifurcation of the splay phase when cos α = 0, namely at α = π/2 independent of β. There can also be bifurcation from splay phase to rotating block where meaning there is a bifurcation from the splay phase to rotating blocks (0, π, θ, θ + π) when cos β = 0. The dynamics of the points with isotropy (S 1 × S 1 ) × s Z 2 (rotating blocks (0, π, θ, θ + π)) is determined bẏ θ = 4r cos θ sin(θ − β) meaning this is degenerate for β = π/2, independent of α! Note that the bifurcations for the rotating blocks only becomes nondegenerate on addition of the fourth harmonic-only then can we get nontrivial solutions with this symmetry. There are two possible two-cluster state: for N = 4 and p = 1 so that P = 1/4 and setting q = −1, (Equation 23) can be expressed as the cubic equation  with discriminant 16 (cos (β)) 2 r 2 − 4 (cos (α)) 2 = 0.
Hence the bifurcations of these states are at The remainder of the curves on Figure 4 are obtained by numerical path following of local bifurcations on the two dimensional fixed point subspace of S 2 in the reduced system (Equation 11) using XPPAUT [28]. We remark there may well be more local bifurcations that involve fully symmetry broken (fourcluster) states, as well as global bifurcations that are not shown on this diagram.

DYNAMICS FOR ODD AND EVEN COUPLING FUNCTIONS
In the case that the coupling function g in Equation (1) is even or odd, the system has additional structure. Write g = g + + g − where are the even and odd parts respectively. For the coupling functions of the form (Equation 4) we have g + (ϕ) = −q cos(ϕ) sin(α) − r cos(2ϕ) sin(β), g − (ϕ) = q sin(ϕ) cos(α) + r sin(2ϕ) cos(β).
The remaining results in this section are independent of N and apply to general even or odd g. -see there for more details of the dynamics in this case. Note that the case β = π/2 and α = π/2 (or 3π/2) correspond to an even coupling function.
If g is odd (i.e., g + (ϕ) ≡ 0) then the system for phase differences is a gradient system [10,29]. More precisely, if we write G ′ (ϕ) = g(ϕ) which is single-valued if g is odd and periodic, we can define a potential then Equation (1) can be written in the form This can be seen as a gradient system for the phase differences (Equation 11) on identifying T N−1 as the set θ ∈ T N such that Hence, the only have attractors that are local minima of the potential V. In general, for ω = 0 one can write a potential for the full system on R N but note that this is multivalued on T N . The case where the coupling function g is even (g − ≡ 0) is dynamically more interesting: there are hints that it organizes the emergence of chaotic dynamics for nearby parameter values [22]. (1) is divergence free. Moreover, the reduced system (Equation 11) is also divergence free.

Lemma 3. If g is even then the flow of Equation
Proof. To see this, consider Frontiers in Applied Mathematics and Statistics | www.frontiersin.org and so Since g is even, this implies that g ′ is odd. Thus, g ′ (θ k − θ j ) = −g ′ (θ j − θ k ) and we have More than this, note that the Jacobian in the direction of the T action is trivial, meaning that the dynamics of the phase differences (Equation 11) transverse to this direction is also divergence free.
Recall that an involution R : that is, R maps any solution to another solution with the direction of the flow is inverted [13].

Lemma 4. If g is even then the system (Equation 1) has a time and parameter-reversal symmetry
For the reduced system (Equation 11) this corresponds to a timereversal symmetry R : Proof. Note that for ω = 0 we have At the same time, which implies that R is a time-reversal symmetry of Equation (1) for the special case ω = 0. This immediately implies that there is a time-reversal symmetry of Equation (11) in general, but only a time and parameter reversal symmetry in general for Equation (1).
Note that given any symmetry g ∈ Ŵ and a time-reversal symmetry R, the composition of the two R • g is also a timereversal symmetry. In fact, the set of time-reversal symmetries on T N−1 is in one-to-one correspondence with the set of symmetries Ŵ [13] and can be generated in this way: given any time-reversal symmetry R, the set of time-reversal symmetries is RŴ. A consequence is that if θ ∈ T N has isotropy H ⊂ S N × T then so does R(θ ).
If R is a time reversal-symmetry, the fixed point subspace Fix(R) = θ ∈ T N | Rθ = θ is not necessarily dynamically invariant, but it is of interest as it yields a natural way to find families of periodic solutions and homoclinic and heteroclinic orbits [13]. In particular, the reversible Lyapunov Centre Theorem asserts that, subject to nonresonance conditions, near an equilibrium in Fix(R) with pairs of purely imaginary eigenvalues there is a two parameter family of reversible periodic orbits, even if there are additional zero eigenvalues [30,31].
There is a particular time-reversal symmetry R : T N−1 → T N−1 that maps the canonical invariant region C to itself. More precisely, let us define on the canonical invariant region and note that this is a reversing symmetry corresponding to reversing the order of the components and composing with R. We characterize all points in C that are mapped to themselves by a reversing symmetry. Thus, we define RC: = θ ∈ C τ q θ = Rθ for some q = 0, . . . , N − 1 (38) These are the fixed points of time-reversal symmetries within C. Let us define and note that On ∂C we have for θ ∈ {0, π} N ∩ C ⊂ RC. Moreover, the splay state (Equation 17) is contained in RC. We now turn to the timereversal symmetries for N = 3 and N = 4 oscillators. In each case we rescale time by N for simplicity.

Dynamics for Even Coupling Functions: N = 3
A system for N = 3 oscillators has and additional degeneracy: there is a constant of motion with G ′ = g which generalizes a specific case of the Watanabe-Strogatz constant of motion [21,32]. Hence, the dynamics of the system is effectively one-dimensional. For the relative fixed point sets of the time-reversal symmetry, we have Q 3,1 = (0, θ 2 , π + θ 2 /2) θ 2 ∈ (0, π) (42b) which define three line segments that meet at the splay state. The configurations of the oscillators are depicted in Figure 5. Figure 6 shows RC = Q N,0 ∪ Q N,1 ∪ Q N,2 superimposed on the dynamically invariant curves determined by the level sets of Equation (41) for an even coupling function with four Fourier modes: The points with time-reversal symmetry RC intersect the boundary of the canonical invariant region at (0, 0, 0), (0, 0, π), (0, π, π).
To find the equilibria in RC it suffices to evaluate the equilibria in Q 3,0 since Q 3,1 , Q 3,2 are the images of Q 3,0 under the action of Z 3 . Write ψ k = θ k − θ 1 and F ψ k = F k − F 1 . On Q 3,0 we have the vector field has eigenvalues where indices are taken modulo N = 3. The eigenvalues thus always take either real or pure imaginary values. Figure 6 shows examples where it appears that the only equilibria in C lie in RC. Equilibria with a pair of imaginary eigenvalues give rise to "islands" of neutrally stable periodic orbits. These are organized by the stable and unstable manifolds of the saddle points on RC which form heteroclinic connections and yield the topological obstructions bounding the islands. The equilibria should relate directly to the constant of motion where the level sets are not smooth manifolds.

Dynamics for Even Coupling Functions: N = 4
For N = 4 we calculate The configuration of phases are depicted in Figure 7. Note that τ Q 4,0 = Q 4,2 and τ Q 4,1 = Q 4,3 . The sets Q 4,1 , Q 4,3 are subsets of planes that intersect the edges of the canonical invariant in the cluster states with S 2 × S 2 isotropy and lines on the faces corresponding to cluster states with S 2 isotropy. The planes intersect in the interior of the canonical invariant region in the set of points with Z 2 isotropy 1 and subdivide the CIR into four connected components. The sets Q 4,0 , Q 4,2 are line segments that intersect ∂C in (0, π, π, π), (0, 0, π, 0), (0, 0, 0, π), and (0, π, 0, 0). We have 3 j = 0 Q 4,j = splay . Figure 8 shows a projection of the canonical invariant region and RC for coupling function (Equation 43) with four nontrivial Fourier modes. Note that while RC is not necessarily invariant as fixed point subspaces of a time-reversal symmetry, it may contain invariant subsets, such as points of nontrivial isotropy.
Similarly calculating the eigenvalues, we find and the eigenvector for λ = 0 is given by Since g is even, there are solutions to the fixed point (Equation 51) that are independent of g. The arguments being equal, ψ ⋆ 2 −ψ ⋆ 3 = ψ ⋆ 2 + ψ ⋆ 3 + 2qπ, q ∈ Z, implies ψ ⋆ 3 = qπ, q ∈ Z and arbitrary ψ ⋆ 2 . If the sign of the arguments is opposite, we have The condition ψ ⋆ 3 = π defines a line of equilibria given by L − = (0, ϕ, π, π + ϕ) ϕ ∈ (0, π) which are points of Z 2 isotropy in C. Note that the eigenvector v 0 of the trivial eigenvalue, Equation (53), points along L − . Moreover, the remaining two eigenvalues evaluate to If some harmonics of the coupling function g vanish they take either real or imaginary values for all points in L − : if the coupling function is π-periodic, g(φ + π) = g(φ), L − is an continuum of saddles since λ = ±2g ′ (ϕ).
If g(φ + π) = −g(φ) we obtain implying that neutrally stable periodic orbits arise close to L − . Similarly, in the second case where the sign of the arguments is opposite, consider ψ ⋆ 2 = 0 which yields a continuum of equilibria which are points of S 2 × S 2 isotropy on ∂C. Again, the eigenvector v 0 corresponding to the neutrally stable direction points along L + . The remaining eigenvalues are given by yielding a family of saddles. The stable and unstable manifolds lie in boundary faces of C with S 2 isotropy. They can intersect RC again in Q 4,1 which intersects the boundary faces adjacent to L + in a line. While the equilibria studied above exist for any choice of even coupling function g, there may be further equilibria that organize the dynamics. First, the Equations (47), (51) may have more equilibria that depend on the coupling function chosen. As these lie in RC they have symmetric stability properties. Second, the equilibria that are not contained in RC appear to play an important role and yield a mechanism that could give rise to complex heteroclinic networks in the system. For example, assume that θ ⋆ ∈ C RC is a hyperbolic equilibrium with trivial isotropy. Since the vector field is divergence free, the equilibrium has to be a saddle. Suppose that θ ⋆ has onedimensional stable and two-dimensional unstable manifold. The action of symmetries and time-reversal implies that there are N symmetric copies of θ ⋆ . At the same time, there are N copies with inverse stability properties that are given by the group orbit of R(θ ⋆ ). An intersection of the unstable manifold of θ ⋆ with Q 4,1 or Q 4,3 yields a heteroclinic connection to an equilibrium related by the time-reversal symmetry. Such a heteroclinic connection may have trivial isotropy (given the way Q 4,1 and Q 4,3 subdivide the CIR into connected components) or nontrivial isotropy, depending on where this intersection takes place. In this way, a network between the (symmetrically and reversing symmetry related) copies of θ ⋆ may be formed.

Even Coupling and Integrability
As noted in Equation (41) for N = 3 and arbitrary even coupling function, there is N − 2 = 1 constant of the motion, in agreement with the case for Kuramoto coupling [15]. It is natural to consider whether there are N − 2 independent constants of the motion for N = 4 and arbitrary even coupling function. If this was the case, the only invariant sets would be two frequency quasiperiodicity. Although Watanabe and Strogatz [15] and the presence of time-reversal symmetries might hints this might be the case, we present evidence there is at most one integral of the motion for N = 4.
In particular, for N = 4 and the states with S 2 isotropy on the boundary of C, Figure 9 demonstrates that for the case N = 4 and even coupling functions with up to four harmonics of the form (Equation 43) are not fully integrable. In cases (c), (d) of the Figure, we argue there are no nontrivial smooth integrals of the motion within the S 2 subspace, owing to the appearance of sink/source pairs. We only find such evidence of nonintegrable behavior if c 3 and c 4 are non both zero: if c 3 and c 4 are both zero (Equation 43) (only two harmonics) then there are apparently no equilibria within S 2 that are not fixed by a time reversal symmetry and the phase portraits are consistent with the presence of a smooth constant of the motion.
It is still possible that there is one integral of the motion for the case N = 4 and general even coupling, as long as this integral is zero on the boundary of the canonical invariant region, and the trajectory shown in Figure 8 is consistent with this.

Open Questions
We finish by highlighting some open questions raised in this paper. Firstly, we note that addition of the second harmonic does unfold many degeneracies, though not all. We note that there may be some that are still present for (Equation 1) even for coupling functions with an infinite number of harmonics present, and these will disappear on considering more general coupling that is not of "pairwise" type (see [8]). There seems to be a subtle interplay of the order N, the symmetries of g and the number of harmonics L that determine the dynamics.
Extending a bifurcation analysis to higher N will be harder than extending to higher L, but the latter seems worthwhile especially if it is possible to determine the location of global bifurcations. We do extend the bifurcation analysis here in part analytically, using the special form of the coupling function to write the bifurcations of two-cluster in polynomial form. Potentially this can be extended to more general cluster states, though at the expense of ever more complex expressions. As an alternative approach, rather than doing a bifurcation analysis for fixed L, one could conversely ask whether some harmonics are crucial for specific bifurcations to occur.
The dynamics for the case of an even coupling function shows remarkable subtlety. Numerical experiments show a variety of periodic orbits which can be quite complicated in examples of even coupling functions where a higher harmonic dominates; cf. Figure 8. Although we demonstrate in the previous section that one cannot expect there to be fully integrability with N − 2 independent constants of the motion [15,Theorem,p. 249], there may well be one constant of the motion in the general even coupling case which we have not yet been able to determine.
For even coupling we demonstrate for N = 3, 4 that the set of phases with time-reversal symmetry within the canonical invariant region is a nontrivial (typically noninvariant) set that organizes the dynamics. Complicated dynamics, such as chimera states [23] and chaotic dynamics [22] (and as a consequence also the chaotic weak chimeras constructed in Bick and Ashwin [33]), typically arise in systems for small perturbations away from a system with time-reversal symmetry. We anticipate that understanding these heteroclinic networks in greater detail will shed light on how complicated arise in networks of coupled phase oscillators. Naturally, the more equilibria there are, the more intricate the heteroclinic network can be. Thus, deriving conditions for the existence of equilibria depending on the coupling function seems like a first step to understand the dynamics for even coupling functions g.