Heteroclinic Dynamics of Localized Frequency Synchrony: Heteroclinic Cycles for Small Populations

Many real-world systems can be modeled as networks of interacting oscillatory units. Collective dynamics that are of functional relevance for the oscillator network, such as switching between metastable states, arise through the interplay of network structure and interaction. Here, we give results for small networks on the existence of heteroclinic cycles between dynamically invariant sets on which the oscillators show localized frequency synchrony. Trajectories near these heteroclinic cycles will exhibit sequential switching of localized frequency synchrony: a population oscillators in the network will oscillate faster (or slower) than others and which population has this property changes sequentially over time. Since we give explicit conditions on the system parameters for such dynamics to arise, our results give insights into how network structure and interactions (which include higher-order interactions between oscillators) facilitate heteroclinic switching between localized frequency synchrony.


Introduction
Networks of interacting oscillatory units can give rise to dynamics where the system appears to be in one metastable state before "switching" to another in a rapid transition. Such dynamics are in particular believed to be of functional relevance for neuronal networks where one observes sequential switching between patterns involving localized activity or synchrony [1,2,3]. One approach is to capture these dynamics on a macroscopic scale: one assigns each pattern an activity variable whose dynamics are then described by kinetic equations [4]. The resulting equations are of generalized Lotka-Volterra type which support stable heteroclinic cycles, that is, cycles of hyperbolic equilibria which are connected by heteroclinic trajectories. The dynamics near such heteroclinic cycles now resemble sequential switching dynamics of activity patterns. Indeed, heteroclinic cycles and networks have been long studied in their own right; see [5] for a recent review.
However, such a qualitative approach fails to capture the dynamics on the level of single, nonlinearly interacting oscillators. In particular, is does not necessarily illuminate what ingredients of network topology and the interactions between oscillators [6] facilitate switching dynamics. If one assumes weak coupling, phase reduction provide a powerful tool to describe the dynamics of an oscillator network; in this reduction, each oscillator is represented by a single phase variable on the torus T := R/2πZ and the Date: October 17, 2018. 1 dynamics of the phases are described by a phase oscillator network. Simple networks of globally and identically coupled phase identical oscillators support heteroclinic cycles and networks [7]. The equilibria involved in these cycles are phase-locking patterns with oscillators in different clusters which have a constant phase difference. The symmetry properties of these networks, however, imply that all oscillators rotate with the same speed (frequency) on average-the network is globally frequency synchronized. In a neural network, this corresponds to all neurons to fire at the same average rate while the exact timing of firing changes.
By contrast, even networks of identical phase oscillators that are organized into different populations can give rise to dynamics where frequency synchrony is local to a population rather than global across the whole network. In other words, the interactions in a network of identical oscillators cause some units to evolve faster (or slower) than others. Dynamically invariant sets with this property relate to "chimeras" [8,9,10] who have-as patterns withe localized frequency synchrony-been hypothesized to play a functional role in the context of neuroscience [11,3,12]. From a mathematical point of view, the notion of a weak chimera [13,14,15] formalizes the definition of a dynamically invariant set with localized frequency synchrony for finite networks of identical phase oscillators.
Here we prove the existence of robust heteroclinic cycles between invariant sets with localized frequency synchrony in small phase oscillator networks with higher-order interactions. In contrast to attracting sets with localized frequency synchrony, the dynamics here induce sequential switching dynamics: which population of oscillators oscillates at a faster (or slower) rate will change over time. These results are of interest from several distinct perspectives. First, they illuminate how the interplay of network structure and functional interactions between units gives rise to heteroclinic dynamics in phase oscillator networks: we explicitly relate the network coupling parameters to the existence of heteroclinic cycles. Second, the results highlight how higher-order network interactions shape the (global) network dynamics; apart from higher harmonics in the phase interaction function, the higherorder interactions also include nonadditive interactions between oscillator phases which arise naturally in phase reductions of generically coupled identical oscillators [16] or other resonant interactions [17]. Here, the interplay of higher-order interactions and nontrivial network topology and induces dynamics beyond (full) synchrony. Third, our results provide new examples of heteroclinic cycles in network dynamical systems relevant for applications. We highlight how these examples are distinct from situations previously considered in the literature.
This work is organized as follows. In this paper, we build on results in a recent brief communication [18] to prove the existence of robust heteroclinic cycles between localized frequency synchrony; in a companion paper [19] we give a detailed discussion of the stability of such heteroclinic cycles (which may be embedded into larger heteroclinic structures). The remainder of this paper is organized as follows. In Section 2 we review some preliminaries on heteroclinic cycles and phase oscillator networks. In Section 3 we show existence of a heteroclinic cycle between localized patterns of frequency synchrony in networks consisting of three populations of two oscillators. In Section 4 we consider networks which consist of three populations of three oscillators and show the existence of a heteroclinic cycle of localized frequency synchrony; here, there are continua of saddle connections in two-dimensional invariant subspaces. Finally, in Section 5, we give some numerical evidence that these phenomena persist in networks with more generic interactions before giving some concluding remarks.

Preliminaries
2.1. Heteroclinic cycles. Let M be a smooth d-dimensional manifold and let X be a smooth vector field on M. For a hyperbolic equilibrium ξ ∈ M let W s (ξ) and W u (ξ) denote its stable and unstable manifold, respectively. Definition 2.1. A heteroclinic cycle C consists of a finite number of hyperbolic equilibria ξ k ∈ M, k = 1, . . . , Q, together with heteroclinic trajectories where indices are taken modulo q.
For simplicity, we write C = (ξ 1 , . . . , ξ Q ). If M is a quotient of a higherdimensional manifold and C is a heteroclinic cycle in M, we also call the lift of C a heteroclinic cycle.
While heteroclinic cycles are in general a nongeneric phenomenon, they can be robust in dynamical systems with symmetry. Let Γ be a finite group which acts on M. For a subgroup H ⊂ Γ define the set Fix(H) = { x ∈ M | γx = x ∀γ ∈ H } of points fixed under H; any Fix(H) is invariant under the flow generated by X. For x ∈ M let Γx = { γx | γ ∈ Γ } denote its group orbit and Σ(x) = { γ ∈ Γ | γx = x } its isotropy subgroup. Now assume that the smooth vector field X is Γ-equivariant vector field on M, that is, the action of the group commutes with X. Now let C = (ξ 1 , . . . , ξ Q ) be a heteroclinic cycle with the following properties. For an isotropy subgroup Σ q ⊂ Γ write P q = Fix(Σ q ). Now suppose that there exist Σ q (and thus P q ) such that ξ q , ξ q+1 ∈ P q , ξ q+1 is a sink in P q , and [ξ q → ξ q+1 ] ⊂ P q . Then C is robust with respect to Γ-equivariant perturbations of X, that is, Γ-equivariant vector fields close to X will have a heteroclinic cycle close to C; see [20] for details.
2.1.1. Dissipative heteroclinic cycles. Trajectories close to a heteroclinic network can show switching dynamics: qualitatively speaking, the trajectory will spend time close to one saddle ξ q before a rapid transition to ξ q+1 . This is in particular the case when the heteroclinic cycle is attracting (in some sense); see for example [5] for a more elaborate discussion.
Here we consider a criterion where we expect attraction in some sense based on the local attraction properties at a hyperbolic equilibrium ξ. Let λ (j) denote the eigenvalues of the linearization of X at ξ ordered such that Re λ (1) ≤ · · · ≤ Re λ (l) < 0 < Re λ (l+1) ≤ · · · ≤ Re λ (d) .
Intuitively speaking, a heteroclinic cycle is dissipative if there is more contraction of phase space than expansion close to the saddle points. Obviously, a heteroclinic cycle is dissipative if all its equilibria are dissipative. The following result shows that, subject to suitable additional assumptions, we may expect a dissipative heteroclinic cycle to be asymptotically stable.
Remark 2.3. Suppose that Γ is finite and let C be a dissipative heteroclinic cycle in R d such that In other words, the entire unstable manifold of one saddle is contained in the stable manifold of the next saddle. Then the results in [23] imply that C is asymptotically stable.
Here we restrict ourselves to show the existence of dissipative heteroclinic cycles; we address the problem of stability explicitly in the companion paper [19]. Indeed, if there are unstable manifolds of dimension larger than one, there may be continua of heteroclinic trajectories. In such a case, the question about stability is more challenging as discussed in [24], in particular because the condition in Remark 2.3 does not allow any set of points (however small) on the unstable manifold of one saddle to lie outside of the stable manifold of the next saddle.
We recall some definitions given in [24], adapted to our setting. Suppose that C is a heteroclinic cycle. We say that there is a (directed) edge between This defines a directed graph G(C). LetG(C) := G(C)/Γ denote the quotient obtained by identifying vertices and connections on the same group orbits. The graphG(C) is cyclic if each vertex has unique edges entering and leaving it.
Definition 2.4. For a heteroclinic cycle C define the associated heteroclinic chain IfG(C) is cyclic then H(C) is a cyclic heteroclinic chain.
In contrast to Definition 2.1, the heteroclinic chain associated to a heteroclinic cycle now contain all heteroclinic trajectories which connect individual equilibria. Note that heteroclinic chains do not need to be closed in M: some part of W u (ξ q ) for some q may lie outside of the heteroclinic chain.
2.2. Phase oscillator networks with nonpairwise interactions. Consider M populations of N phase oscillators where θ σ,k ∈ T denotes the phase of oscillator k in population σ. Hence, the state of the oscillator network is determined by θ = (θ 1 , . . . , θ M ) ∈ T N M where θ σ = (θ σ,1 , . . . , θ σ,N ) ∈ T N is the state of population σ. Let g 2 , g 4 : T → R be 2π-periodic functions and Now consider the phase oscillator network where the phases of individual oscillators evolve according tȯ where ω is the intrinsic frequency of each oscillator 1 . For these network dynamics, the phase interactions within populations are determined by the coupling (or phase interaction) function g 2 whose arguments differences of oscillator pairs. By contrast, the interactions between populations, given by (1), are mediated by the nonpairwise interaction function g 4 whose argument is a linear combination of four of the oscillators' phases. The parameter K − > 0 determines the coupling strength to the previous population whereas K + > 0 determines the coupling strength to the previous population. Here we assume K := K − = K + > 0 for simplicity. For g 4 = cos, the equations (2) approximate the dynamics of a phase oscillator networks with mean-field mediated bifurcation parameters up to rescaling of time as outlined in [18].
corresponds to phases being in full phase synchrony and denotes a splay phase configuration-typically we call any element of the group orbit S N D a splay phase. For a network of interacting populations, we use the shorthand notation to indicate that population σ is fully phase synchronized or in splay phase. Consequently, S · · · S (M times) is the set of cluster states where all populations are fully phase synchronized and D · · · D the set where all populations are in splay phase.
The network interactions in (2), which include nonpairwise coupling, induce symmetries. More precisely, the equations (2) are (S N × T) M Z Mequivariant. Each copy of T acts by shifting all oscillator phases of one population by a common constant while S N permutes its oscillators. The action of Z M permutes the populations cyclically. These actions do not necessarily commute.
The symmetries yield invariant subspaces on T M N for the dynamics given by (2). In particular, the S N permutational symmetries within each population imply that the sets (5) are invariant [25]. Moreover, any set of the form θ 1 · · · θ M with θ k ∈ {S, D} is an equilibrium relative to the continuous T M symmetry, that is, the corresponding ψ 1 · · · ψ M is an equilibrium in the reduced dynamics.

2.2.2.
Frequencies and localized frequency synchrony. Suppose that M > 1 and let θ : [0, ∞) → T M N be a solution of (2) with initial condition θ(0) = θ 0 . Whileθ σ,k (t) is the instantaneous angular frequency of oscillator (σ, k), define the asymptotic average angular frequency of oscillator (σ, k) by Here we assume that these limit exists for all oscillators but this notion can be generalized to frequency intervals; see also [14,15].
localized frequency synchrony if for any θ 0 ∈ A we have Ω σ,k = Ω σ and there exist indices σ = τ such that Remark 2.6. Note that a chain-recurrent set A with localized frequency synchrony is a weak chimera as defined in [13].

Heteroclinic Cycles for Two Oscillators per Population
In this section, we show the existence of robust heteroclinic cycles for networks of M = 3 populations of N = 2 oscillators. Let g 2 , g 4 : T → R denote the 2π-periodic coupling functions which govern the interactions within and between populations as above. With the network dynamics (2) can be written aṡ By reducing the T M symmetry, we obtain the dynamics of the phasedifferences aṡ The phase space of (9) is organized by invariant subspaces as sketched in Figure 1. For completeness, we characterize SSS, DDD before we focus on sets with localized frequency synchrony.
In the reduced dynamics (10), both SSS and DDD are equilibria. The equilibrium SSS = (0, 0, 0) lies in the intersection of the invariant subspaces ψSS, SψS, and SSψ. On these subspaces, the dynamics are given bẏ Thus, the linear stability of SSS is determined by the triple eigenvalue λ SSS = −2g 2 (0) (12) which correspond to a perturbation separating the phases of one population. Similarly, DDD = (π, π, π) lies in the intersection of the invariant subspaces ψDD, DψD, and DψD. On these invariant subspaces, the dynamics are determined by (13)ψ = 2ĝ 2 (ψ), as well. Linearizing at ψ = π yields (14) λ DDD = −2g 2 (π) which determines the stability of DDD. In the full system (9), there are three additional zero eigenvalues for eigenvectors along the group orbit of the phase-shift symmetry. Note that the linear stability of SSS, DDD are fully determined by the pairwise coupling g 2 within populations.
3.1. Saddle invariant sets with localized frequency synchrony. For the remainder of this section, we will consider (8) with interaction given by the coupling functions These interactions include higher harmonics; in particular, for r > 0 the calculations above imply that both SSS and DDD are linearly stable. Moreover, for α := α 2 = α 4 − π 2 we obtain the same parametrization as in [18,19].  (13), on broken lines by (16), and on solid lines by (17). Panel (b) shows a heteroclinic cycle between these equilibria.
In the following we estimate the asymptotic average frequencies of DSS, DDS, SDS, SDD, SSD, DSD, and DSS. Note that it suffices to consider DSS, DDS since the latter four are their images under the Z M action which permutes populations.
In this section, we show that there are dissipative robust heteroclinic cycles for networks of M = 3 populations of N = 3 oscillators with coupling (24). Indeed, we proceed as before and derive conditions for which there heteroclinic source-sink connections on invariant subspaces forced by symmetry. The resulting heteroclinic network will be robust. Thus, for the remainder of the section we set (α 2 , α 4 ) = ( π 2 , π) rather writing down the conditions in full generality.
A stability analysis of the phase configurations SSS and DDD can be done as in the previous section. Moreover, we restrict ourselves to DSS and DDS because of symmetry. 4.1. Local dynamics. We first establish conditions for DSS, DDS to be invariant sets with localized frequency synchrony which are suitable saddles in the reduced system. Proof. The proof is essentially the same as for Lemma 3.1. Frequency synchrony with populations is given by Lemma 2.7. For K = 0 we have Ω 1 (θ 0 ) = ω + 2g 2 (0) = ω + 2 for θ 0 ∈ Sψ 2 ψ 3 and Ω 2 (θ 0 ) = ω + g 2 (2π/3) + g 2 (4π/3) = ω − 1 for θ 0 ∈ ψ 1 Dψ 3 . This implies that for K ≥ 0 and coupling (24)   Similarly, for the linearization of the vector field at DDS we have which govern the linear stability of the phase configuration of the first, second, and third population, respectively.
Thus, we have hyperbolic saddles with two-dimensional unstable manifold if 0 < r and which is equivalent to (CλN3), as asserted.
Note that expansion at DSS is determined by the double real eigenvalue λ DSS 2 forced by symmetry. The eigenvalues of the linearization (25) and (26) also give insight into the local bifurcations as parameters are varied. For example, DDS undergoes a Hopf bifurcation as the parameter r goes through zero.

Global dynamics.
In the previous section we established the existence of suitable saddle invariant sets. To obtain a heteroclinic cycle, these invariant sets have to be joined by a heteroclinic trajectory.
Both DSS and DDS lie in the two-dimensional invariant subspace DψS and DDS while SDS lie in the two-dimensional invariant subspace ψDS. Because of the permutational symmetry within populations, the dynamics on both DψS and ψDS are S 3 equivariant. As discussed in the context of Lemma 2.7, this implies that the phase ordering within each population is preserved. Hence it suffices to consider the dynamics on (the closure of) the invariant simplex, commonly referred to as the canonical invariant region, for the phase differences; cf. [25,26]. The phase configuration S = (0, 0), where all oscillators are phase synchronized, lies on the boundary ∂C of C and the splay phase configuration D = 2π 3 , 4π 3 is its centroid as illustrated in Figure 3(a). For a function F : C → R let F C denote the sup norm on C.
The vector field for the dynamics of (2) with coupling (24) and α 2 = π 2 , α 4 = π can now be evaluated explicitly. Write The dynamics on DψS are given by and the dynamics on ψDS are given by The first term describes interaction within each population given by the first harmonic only, the second term is the intra-population interaction through higher harmonics, and the third term are interactions arising through the coupling between populations.  Proof. First note that the dynamics of the uncoupled system, K = 0, without higher harmonics, r = 0, is integrable (on both DψS and ψDS) since the quantity is preserved; see also [26]. (Some level sets of V are depicted in Figure 3(a).) That is, we haveV := grad V,ψ = grad V, X 0 = 0.
Note that V > 0 on C, it vanishes on its boundary ∂C, and takes a unique maximum V (D) at D. We will first consider the dynamics on DψS and show that there is a source-sink heteroclinic trajectory [DSS → DDS]. By assumption (CλN3), we have that DSS is a source in DψS with W u (DSS) ⊂ DψS. In the following, we derive conditions on K and r such that V is a (Lyapunov-like) potential function on C which guarantee that V is strictly increasing along trajectories in C. Thus, any trajectory in C converges to D ∈ C which yields a heteroclinic trajectory [DSS → DDS]. Now consider nontrivial coupling between populations, K > 0, while higher harmonics are absent, r = 0. Using trigonometric identities we have The potential V increases along trajectories ifV > 0 on C D. Rearranging terms,V > 0 is equivalent to for ψ ∈ C D. The function Q has a unique maximum on C at ψ = D with Q(D) = 1. Thus,V > 0, and trajectories in C approach D asymptotically. We now show that this property persists for r > 0 sufficiently small. We haveV = K grad V, X K + r grad V, X r and by the calculations above, we know that grad V, X K only vanishes for ψ ∈ ∂C ∪ D. The conditionV > 0 on C D is equivalent to Note that any singularity of R is removable since grad V, X r also vanishes on C D at the same order. Hence, R(ψ) is a bounded function on C and evaluating minima and maxima on C yields R C < 15; cf. Figure 3(b). This implies thatV > 0 on C D if 15r < K.
which yields a heteroclinic connection between the source S and the sink D on DψS. The proof that there is a robust heteroclinic trajectory [DDS → SDS] in ψDS is analogous. We have W u (DDS) ⊂ ψDS and a heteroclinic trajectory [DDS → SDS] is a trajectory connecting the source D and sink S on ψDS. Thus, it suffices to show thatV < 0 on C D. Note that the vector fields (28) and (29) only differ by the sign of K. We proceed as above and in the last step the conditionV < 0 is equivalent to Again, R C < 15 implies that if 15r < K we have a robust heteroclinic connection between the source D and the sink S on ψDS.
In fact, Lemma 4.3 implies the existence of a heteroclinic cycle  Proof. In order to get dissipativity of the cycle, we need to control the product of the saddle values. More precisely, the cycle is dissipative if It is straightforward to verify that this condition is equivalent to (CνN3).
Note that the conditions (CΩN3), (CλN3), (CψN3), and (CνN3) as given in Lemmas 4.1-4.4 are satisfied simultaneously for (35) 0 Moreover, the heteroclinic trajectories are source-sink connections in invariant subspaces forced by symmetry. This proves the following result. To formalize this observation, we adapt some terminology that was recently introduced in [27].
Definition 4.6. Let C = (ξ 1 , . . . , ξ Q ) be a heteroclinic cycle and H(C) its associated heteroclinic chain. An equilibrium is of measure zero (with respect to any Riemannian measure on W u (ξ q )), (iii) ξ q is equable in H(C) if dim(C pq ) is equal for all p with C pq = ∅. The heteroclinic chain H(C) is complete, almost complete, or equable if all its equilibria are complete, almost complete, or equable respectively.
Note that completeness relates to clean heteroclinic networks defined in [28]. With these notions we obtain the following statement.
Theorem 4.7. For parameters satisfying (35) and r sufficiently small, the heteroclinic chain H(C 3 ) is cyclic, equable, and almost complete but not complete. The closure of H(C 3 ) is complete, but neither cyclic nor equable.
Proof. Note that the heteroclinic chain H(C 3 ) associated with C 3 is cyclicbecause W u (DSS) ⊂ DψS and W u (DDS) ⊂ ψDS-and thus equable.
First, consider the invariant set DψS; it suffices to consider C as above. The invariant set ∂C consists of points of nontrivial isotropy, i.e., Σ(θ) = {id} for θ ∈ ∂C, and ∂C has zero (Lebesgue) measure in C. Since W u (DSS) ∩ ∂C = ∅ and W s (DDS) ∩ ∂C = ∅, the equilibrium DSS cannot be complete. Parametrize a side of ∂C by χ ∈ [0, 2π] (for example χ = ψ 1 = ψ 2 ) where χ ∈ {0, 2π} corresponds to S ∈ ∂C; cf. lie in the same group orbit) with χ ≈ 2 arctan(3K). These are of saddle type, attracting within ∂C and transversely repelling (within DψS); these are shown in Figure 3(c). Therefore, is two-dimensional and of full measure in W u (DSS). This implies that DSS is almost complete. Second, by evaluating the dynamics on ψDS one can show by a similar argument that there are three equilibria ξ ψDS on ∂C. These are attracting transversely to ∂C and repelling within ∂C. Consequently, and DDS is almost complete. However, DDS is not complete since W s (ξ ψDS ) ⊂ W u (DDS) is not contained in W s (SDS).
The heteroclinic chain H(C cl ) associated with C cl contains H(C 3 ) but is not cyclic nor equable. Indeed, the graphG(C cl ) for H(C cl ) contains the noncyclic subgraph depicted in Figure 4. This completes the proof of the assertion.

Dynamics of Networks with Noise and Broken Symmetry
The heteroclinic cycles lead to switching between localized frequency synchrony which are observed in numerical simulations. First, define the Kuramoto order parameter of population σ as Z σ = 1 N N j=1 exp(iθ σ,k ) and let R σ = |Z σ |. In particular, R σ ∈ [0, 1] encodes the level of synchrony in each population, that is, R σ = 1 iff θ σ ∈ S and R σ = 0 if θ σ ∈ S. Write (2) aṡ θ σ,k = X σ,k (θ), setting ω = −(N − 1)g 2 (0) without loss of generality such that oscillators appear stationary if phase synchronized in the absence of interpopulation coupling. Let W σ,k be independent Wiener processes with mean zero and variance one. Since attracting heteroclinic cycles show exponential slowing down of transition times between subsequent saddles, we solve the stochastic differential equation for η > 0 using XPP [29]. As shown in Figure 5, the dynamics exhibit switching between localized frequency synchrony: populations sequentially accelerate and decelerate as the populations synchronize and desynchronize. The transition times scale with the noise strength as expected [30]. From the point of view of a phase reduction of a general network of nonlinear oscillators, the interaction terms in the phase oscillator network (2) are nongeneric. Indeed, one would expect that nontrivial pairwise interaction terms would not only be present within populations but also between populations. Here, we asses the effect of forced symmetry breaking on the dynamics (2) in numerical simulations. More specifically, define The function Y sym σ,k is S M N equivariant and yields pairwise interactions between oscillators between different populations. Let ι(σ, k) = (σ − 1)N + k be a linear indexing for all oscillators. Define which yields additional pairwise interaction terms. For parameters δ sym , δ asym we now consider the evolution of (37)θ σ,k = X σ,k (θ) + δ sym Y sym σ,k (θ) + δ asym Y asym σ,k (θ) + ηW σ,k . Heteroclinic switching dynamics persist if the phase-shift symmetries are broken. For δ sym > 0, δ asym = 0, the (S N × T) M Z M symmetry of (2) is broken; if Γ = (S M N Z M ) × T then (37) is Γ-equivariant. In other words, rather than having a phase-shift symmetry for each population, the system (37) has a single phase-shift symmetry which acts by adding a constant phase to all M N oscillators. While this breaks the invariant subspace structure that gave rise to the robust heteroclinic cycles, we expect certain normally hyperbolic tori to persist for δ sym > 0 sufficiently small. More precisely, if DSS = (θ 1 , θ 2 , θ 3 ) ∈ T M N θ 1 ∈ D, θ 2 , θ 3 ∈ S, θ 2 = θ 3 ⊂ DSS (and SDS, SSD are its images under the Z M symmetry action) then we expect that invariant tori exist close to DSS, SDS, and SSD. Indeed, solving the system numerically-as shown in Figure 6-indicates that there is in fact a residual attracting heteroclinic network which approaches DSS, SDS, and SSD.
With further forced symmetry breaking, δ sym , δ asym > 0, the phase oscillator network (37) exhibits irregular switching of localized frequency synchrony even in the absence of noise. These potentially chaotic dynamics arise close to the heteroclinic networks C N , N = 2, 3, as shown in Figure 7.

Discussion and Conclusions
Phase oscillator networks with higher-order interactions can give rise to heteroclinic cycles between frequency synchrony; in numerical simulations, these lead to sequential acceleration and deceleration of oscillator populations. Indeed, because of dissipativity, we expect that the attractor of the deterministic system is a subset of the closure of the associated heteroclinic chain. For networks of N = 2 oscillators in each population we calculate the stability of the heteroclinic cycles and their bifurcations explicitly in the companion paper [19]. For N = 3 the unstable manifold of each saddle is (at least) two-dimensional and the assumptions to apply existing stability results [23,24] are not satisfied. We will address this question in future research. Switching between localized frequency synchrony persists for the network (37) as the phase shift symmetries are broken, δ sym = 0.1; the other parameters are as in Figure 5. Due to the attractive coupling between populations, the synchronized populations now synchronize in phase with each other. In other words, trajectories approach DSS, SDS, and SSD cyclically. . Further symmetry breaking leads to deterministic, irregular cycling of localized frequency synchrony for the dynamics of (37). Here δ sym = 0.1 and δ asym = 0.3δ sym while noise is absent, η = 0; all other parameters are as in Figure 5.
Rather than assuming weak coupling between populations, the results presented here rely on the symmetries induced by the nonpairwise higherorder network interaction terms. Our numerical simulations for nearby vector fields where these symmetries were broken indicated the persistence of some residual heteroclinic structure. In this context, it would be desirable to extend the methods of forced symmetry breaking [31,32,33] to understand the bifurcation behavior for nearby network vector fields with generic interactions.
Numerical simulations indicate that switching dynamics between localized frequency synchrony also arises in networks with M = 3 populations of N > 3 phase oscillators [18]. Indeed, the methods used here are likely applicable to such networks as well: without higher harmonics, r = 0, the oscillators are sinusoidally coupled and the phase space T M N is foliated by low-dimensional manifolds [34,35] on which we expect to have a similar potential functions as in the proof of Lemma 4.3. While we cannot expect hyperbolicity in this limit due to the degeneracy in the system, suitable network interaction terms with higher harmonics, combined with an approach similar to [36], could give rigorous results to show the existence of heteroclinic networks.
In summary, heteroclinic switching dynamics between localized frequency synchrony may arise in networks of identical phase oscillators with higherorder interactions. Here we gave rigorous results for small oscillator networks, but we anticipate similar approaches to be viable for larger networks. While the heteroclinic switching observed here are distinct from those discussed in [37], where large networks (N ≥ 1000) of nonidentical oscillators are considered, it would be interesting to relate the two.