Heteroclinic Dynamics of Localized Frequency Synchrony: Stability of Heteroclinic Cycles and Networks

In the first part of this paper, we showed that three coupled populations of identical phase oscillators give rise to heteroclinic cycles between invariant sets where populations show distinct frequencies. Here, we now give explicit stability results for these heteroclinic cycles for populations consisting of two oscillators each. In systems with four coupled phase oscillator populations, different heteroclinic cycles can form a heteroclinic network. While such networks cannot be asymptotically stable, the local attraction properties of each cycle in the network can be quantified by stability indices. We calculate these stability indices in terms of the coupling parameters between oscillator populations. Hence, our results elucidate how oscillator coupling influences sequential transitions along a heteroclinic network where individual oscillator populations switch sequentially between a high and a low frequency regime; such dynamics appear relevant for the functionality of neural oscillators.


Introduction
Interacting populations of identical oscillators are capable of generating global dynamics that exhibit rapid transitions between metastable states where different populations are in different frequency regimes. Such dynamics can be caused by trajectories close to heteroclinic structures between invariant sets where frequency synchrony is local rather than global across all populations [1]. In the first part of this paper [2], we showed the existence of heteroclinic cycles in three coupled small populations.
For such dynamics to be observable over longer timescales, the heteroclinic cycles have to have some stability properties. Apart from asymptotic stability and instability, heteroclinic cycles can display various intermediate forms of nonasymptotic attraction. These range from fragmentary asymptotic stability ("attracting more than nothing") to essential asymptotic stability ("attracting almost everything"). Podvigina and Ashwin [3] introduced a stability index to quantify attraction along trajectories. This stability index is defined for any dynamically invariant set and thus provides a convenient tool to describe the stability of heteroclinic trajectories within a cycle or network 1 . Recently, Garrido-da-Silva and Castro [4] derived explicit expressions for the stability indices for a fairly general class of heteroclinic cycles called quasi-simple. Such expressions are particularly useful to describe the stability of heteroclinic cycles that are part of a networks consisting of more than one cycle.
The main contributions of this paper are explicit stability results for heteroclinic cycles and networks between invariant sets with localized frequency synchrony in terms of the coupling parameter of the oscillator populations.
Here we focus on coupled oscillator populations with two oscillators per population. Due to the existence of invariant subspaces, the heteroclinic cycles are quasi-simple. Consequently, we apply the stability results of Garridoda-Silva and Castro [4] to calculate the stability indices. We first consider three coupled oscillator populations to calculate stability indices for the heteroclinic cycles in [2]. We then show that four coupled oscillator populations support a heteroclinic network which contains two distinct heteroclinic cycles of the type considered before. Their stability properties are then calculated using the tools developed for three populations and we comment on the stability of the whole network. Since our stability conditions explicitly depend on the coupling parameters of the oscillator populations, our results elucidate how the coupling structure of the system shapes the asymptotic dynamical behavior. Moreover, they highlight the utility of the general stability results in [4] for heteroclinic cycles on arbitrary manifolds.
The remainder of this paper is structured as follows. The following section summarizes facts on (robust) heteroclinic cycles, nonasymptotic stability, and coupled populations of phase oscillators. In Section 3 we calculate the stability indices along the heteroclinic cycle in the first part of this paper [2] for a system of three populations. Such cycles are contained in a heteroclinic network for four coupled populations as shown in Section 4, and we calculate their stability properties. We also give some numerical results and comment on the stability of the network as a whole. Finally, we give some concluding remarks in Section 5.

Preliminaries
To set the stage, we review some results about heteroclinic cycles, their stability properties, and coupled populations of phase oscillators. In terms of notation, we will follow the first part of the paper [2].
2.1. Heteroclinic cycles and their stability. Let M be a smooth ddimensional manifold and let X be a smooth vector field on M. Define the usual limit sets α(x), ω(x) for the flow on M generated by X and t → ±∞. For a hyperbolic equilibrium ξ ∈ M we write to denote its stable and unstable manifold, respectively.
Definition 2.1. A heteroclinic cycle C consists of a finite number of hyperbolic equilibria ξ q ∈ M, q = 1, . . . , Q, together with heteroclinic trajectories where indices are taken modulo Q.
A heteroclinic network N is a connected union of two or more distinct heteroclinic cycles.
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. The same goes for a heteroclinic network N.
While heteroclinic cycles are in general a nongeneric phenomenon, they can be robust if all connections are of saddle-sink type in (lower-dimensional) subspaces. Let C = (ξ 1 , . . . , ξ Q ) be a heteroclinic cycle. If there are flowinvariant submanifolds P q ⊂ M such that [ξ q → ξ q+1 ] ⊂ P q is a saddlesink connection, then C is robust with respect to perturbations of X which preserve the invariant sets P q .
Robust heteroclinic cycles may arise for example 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; this is a vector space that 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 a Γ-equivariant vector field on M, that is, the action of the group commutes with X. Any heteroclinic cycle with P q = Fix(Σ q ) where Σ q are isotropy subgroups is robust to Γ-equivariant perturbations of X, that is, Γequivariant vector fields close to X will have a heteroclinic cycle close to C; see [5] for more details.
2.1.1. Nonasymptotic stability. Heteroclinic cycles may have intricate nonasymptotic stability properties. We briefly recall some definitions that formalize these.
For ε > 0, write B ε (A) for an ε-neighborhood of a set A ⊂ R d and B(A) for its basin of attraction, i.e., the set of points x ∈ R d with ω(x) ⊂ A. For δ > 0 the δ-local basin of attraction is where Φ t is the flow generated by X. Let denote the Lebesgue measure. 6]). An invariant set A is fragmentarily asymptotically stable (f.a.s.) if (B δ (A)) > 0 for any δ > 0.
Being f.a.s. is not necessarily a very strong form of attraction. A set that is not f.a.s. is usually called completely unstable, see also [6]. Melbourne [7] introduces the stronger notion of essential asymptotic stability, which we quote here in the formulation of Brannath [8].

Definition 2.3 ([8], Definition 1.2). A compact invariant set
A is called essentially asymptotically stable (e.a.s.) if it is asymptotically stable relative to a set Υ ⊂ R d with the property that Podvigina and Ashwin [3] introduced the concept of a local stability index σ(x) ∈ [−∞, +∞] to quantify stability and attraction. It is constant along trajectories, so to characterize stability/attraction of a heteroclinic cycle with one-dimensional connections, it suffices to consider finitely many stability indices. Let σ q denote the stability index along [ξ q−1 → ξ q ]. For our purposes it is enough to note that (under some mild assumptions) a heteroclinic cycle C = (ξ 1 , . . . , ξ Q ) is completely unstable if σ q = −∞ for all q, it is f.a.s. as soon as σ q > −∞ for some q, and it is e.a.s. if and only if σ q > 0 for all q = 1, . . . , Q. See [9, Theorem 3.1] for details.
2.1.2. Stability of quasi-simple heteroclinic cycles. The stability indices can be calculated for specific classes of heteroclinic cycles. Let C = (ξ 1 , . . . , ξ Q ) be a robust heteroclinic cycle on M. As above, denote the flow-invariant sets which contain the heteroclinic connections with P q . Let T q := T ξq M denote the tangent space of M at ξ q . For subspaces V ⊂ W ⊂ T q write W V for the orthogonal complement of V in W . In slight abuse of notation, define P − q := T ξq P q−1 and P + q := T ξq P q to be the tangent spaces of P q−1 and P q in T q , respectively. These are linear subspaces of T q of the same dimension as P q−1 (which contains the incoming saddle connection) and P q (containing the outgoing connection), respectively. Set L q := P − q ∩ P + q .
Remark 2.5. Note that this is a slight generalization of the definition given by Garrido-da-Silva and Castro in [4] to arbitrary manifolds. In particular, the condition in Definition 2.4 implies that dim(P q−1 ) = dim(P q ).
As usual, an eigenvalue of the Jacobian dX(ξ q ) is radial if its associated eigenvector is in L q , contracting if the associated eigenvector is in P − q L q , expanding if the associated eigenvector is in P + q L q , and transverse otherwise. In other words, a cycle is quasi-simple if it has unique expanding and contracting directions at each equilibrium, and thus one-dimensional saddle connections.
The standard way to analyze the stability of heteroclinic cycles is to write down a Poincaré return map with linearized dynamics local to the equilibria as well as globally along the connecting orbits; cf. [10]. For quasi-simple cycles whose global maps are rescaled permutations of the local coordinate axes Garrido-da-Silva and Castro [4] showed how their (asymptotic or nonasymptotic) stability can be calculated solely from the properties of the linearization of the equilibria at the cycle. More precisely, the stability of each equilibrium ξ q along the cycle is encoded in a transition matrix M q and the stability of the cycle is determined by properties of these matrices. We explain this technique in more detail when we apply it in Section 3. Note that this immediately implies that the results in [4] carry over to our definition of a quasi-simple heteroclinic cycle since the stability does not depend on other global properties.
For ease of reference, we recall the stability results from [4, Theorems 3.4, 3.10] in a condensed form. For a heteroclinic cycle C = (ξ 1 , . . . , ξ Q ) with transition matrices M q set M (q) := M q−1 · · · M 1 M Q · · · M q+1 M q . All M (q) have the same eigenvalues. If none of the M q has a negative entry-there are no repelling transverse directions-we have the following result, which is a dichotomy between asymptotic stability and complete instability. (i) If M (1) satisfies |λ max | > 1, then σ q = +∞ for all q = 1, . . . , Q and the cycle C is asymptotically stable. (ii) Otherwise, σ q = −∞ for all q = 1, . . . , Q and C is completely unstable.
If the transition matrices M q contain negative entries-there are transversely repelling directions, for example, if the cycle is part of a networkthen additional criteria have to be satisfied in order for the cycle to possess some form of nonasymptotic stability. For a matrix M let λ max denote the maximal eigenvalue and u max = (u max 1 , . . . , u max d ) the corresponding eigenvector. Define the conditions (cf. [4,Lemma 3 > 0 for all m, n = 1, . . . , d. Generally, stability indices are evaluated as a function of the local stability properties at the equilibrium points [3]; for quasi-simple cycles in arbitrary dimension, Garrido-da-Silva and Castro [4] denote this function by F ind . Later on, we will consider three-dimensional transition matrices and for 0 = β = (β 1 , β 2 , β 3 ) ∈ R 3 , this function reads The following proposition summarizes the second stability result adapted to our setting. Here s is bounded by an expression which depends on the number of rows of the transition matrices (and their products) with at least one negative entry; cf. [4] for details.

2.2.
Coupled populations of phase oscillators. Consider M populations of N phase oscillators where θ σ,k ∈ T := R/2πZ denotes the phase of oscillator k in population σ. Hence, the state of the coupled oscillator populations is determined by θ = (θ 1 , . . . , θ M ) ∈ T M N where θ σ = (θ σ,1 , . . . , θ σ,N ) ∈ T N is the state of population σ. Let S N denote the permutation group of N elements. Suppose that the phase evolution is given by where ω is the intrinsic frequency of each oscillator and the vector field Y is Here, each copy of T acts by shifting all oscillator phases of a given population σ by a common constant while S N permutes the oscillator indices k.
The symmetry implies that certain phase configurations are dynamically invariant. For a single population of N oscillators, the subset 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 interacting oscillator 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. Because of the S M N symmetry, the sets (5) are invariant [11].
To reduce the phase-shift symmetry T M we may rewrite (2) in terms of phase differences ψ σ,k := θ σ,k+1 − θ σ,1 , k = 1, . . . , N − 1. Hence, with ψ σ ∈ T N −1 we also write for example ψ 1 S · · · S (or simply ψS · · · S if the index is obvious) to indicate that all but the first population is phase synchronized. The sets (5) are equilibria relative to T M , that is, they are equilibria for the reduced system in terms of phase differences.

2.2.1.
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 limits exist for all oscillators but this notion can be generalized to frequency intervals; see also [12,13].
localized frequency synchrony if for any θ 0 ∈ A and fixed σ we have Ω σ,k = Ω σ for all k and there exist indices σ = τ such that Remark 2.9. Note that a chain-recurrent set A with localized frequency synchrony is a weak chimera as defined in [14].

Three Coupled Oscillator Populations
Here we derive explicit stability results for the heteroclinic cycles in M = 3 coupled populations of N = 2 phase oscillators (2) considered in the first part [2]; we use the same notation introduced there. Interactions between pairs of oscillators are mediated by the coupling function (8) g(ϑ) = sin(ϑ + α) − r sin(a(ϑ + α)).

With the interaction functioñ
the phase dynamics for coupling strength K > 0 between populations are given byθ σ ∈ {1, 2, 3}, k ∈ {1, 2}. These are the equations of motion 2 considered in the first part [2] with phase shifts parametrized by α := α 2 = α 4 − π 2 . The interactions between populations in (10)-which include nonpairwise coupling-are a special case of (2). More precisely, with Z M := Z/M Z the equations (10) are (S N × T) M Z M -equivariant. 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 phase space of (10) is organized by invariant subspaces and there are relative equilibria DSS, DDS and their images under the Z 3 action.
3.1. Heteroclinic cycles and local stability. The coupled oscillator populations (10) with interaction functions (8), (9) support a robust heteroclinic cycle [2]. Linear stability of DSS, DDS are given by the eigenvalues and λ DDS 1 = 2K sin(α) + 2 cos(α) + 4r cos(2α), (12a) For fixed α ≈ π 2 , the assumptions of Lemma 3.1 define a cone-shaped region in (K, r) parameter space: there is an affine linear function L such that K > K 0 where L(K 0 ) = 0 and r between −L(K) and L(K). For the remainder of this section, we assume that the assumptions of Lemma 3.1 hold.
Proof. It suffices to consider the equilibria DSS and DDS due to the symmetry which permutes populations. We have W u (DSS) ⊂ DψS, W u (DDS) ⊂ ψDS which implies that each saddle has one contracting, expanding, and transverse eigenvalue; there are no radial eigenvalues since DSψ ∩ DψS = DSS and DψS ∩ ψDS = DDS.
Subject to nonresonance conditions, we may linearize the flow around the equilibria; see also [15,Proposition 4.1].
(in the second line we allow any combination of + and −). Then we can linearize the flow at the equilibria in C 2 . For α = π 2 , these conditions reduce to r = 0 and r = ±K.
Proof. According to the C 1 linearization theorem [16] we can linearize the flow if the eigenvalues λ l of the linearization satisfy Re λ l = Re λ j + Re λ k when Re λ j < 0 < Re λ k . Given (11), conditions (13) are just these nonresonance conditions. Plugging in α = π 2 yields the second assertion. 3.2. Cross sections, transition matrices, and stability. Using standard notation, we write for the contracting, expanding, and transverse eigenvalues. Thus e q , c q > 0, q ∈ {DSS, DDS}. The ratios between contraction/transverse stability and expansion are given by for q ∈ {DSS, DDS}; we have a q > 0 by definition and b q > 0 if t q < 0.
3.2.1. Poincaré map and transition matrices. We first consider the linearized flow at the equilibria to calculate the local transition maps. Introduce local coordinates (v, w, z) which correspond to the contracting, expanding, and transverse directions, respectively. After appropriate rescaling, consider the cross sections Hence the time of flight is τ = − log(w)/e q which implies that the local map at ξ q is Considering the invariant subspaces, we see that the global maps are rescaled permutations. More specifically, we have Ignoring v, this yields a map between the incoming 2-dimensional sections of subsequent equilibria h q (w, z) = (B q w bq z, D q w aq ).
Taken together, the Poincaré return map for the linearized dynamics around the heteroclinic cycle (modulo the Z 3 group action) is If we introduce logarithmic coordinates we can write the return map in terms of transition matrices [17,4]. Restrict to the (w, z) coordinates and introduce logarithmic variables η = log(w), ζ = log(z). In the new variables, the maps h q become linear, Note that these transition matrices are the same as the ones for simple cycles in R 4 of type C [3].
The transition matrix for the Poincaré map h is M DDS M DSS . These transition matrices govern the stability of the cycle [4, Theorem 3.4].

Stability of
The stability properties at the saddles are symmetric and stability is governed by the properties of the transition matrix Here we omitted the saddle index q since M DSS = M DDS = M. This is the same transition matrix as for a simple heteroclinic cycle in R 4 of type C − 1 [3]. Lemma 3.4 ([3, Section 4.2.2]). A heteroclinic cycle whose stability is given by the transition matrix M is asymptotically stable if b ≥ 0 (that is, t ≤ 0) and a + b > 1; otherwise it is completely unstable.
In terms of the oscillator coupling parameters we can now show that the heteroclinic cycle loses stability completely in a (degenerate) transverse bifurcation at r = 0 as both transverse eigenvalues pass through zero.
Theorem 3.5. For α = π 2 the heteroclinic cycle C 2 is asymptotically stable if r > 0 and completely unstable if r < 0.
Proof. Substituting the stability properties (11), (12), we obtain Simplifying the expressions b ≥ 0 and a+b > 1 now proves the assertion.
These results can now be used to show that the heteroclinic cycle loses stability completely as one of the transverse eigenvalues becomes positive. Proof. First, observe that there are relations between the eigenvalues (11), (12) of the linearization at the saddle points. Set S = 2K sin(α). We have e 1 = S + t 2 , c 1 = S − t 2 , e 2 = S + t 1 , c 2 = S − t 1 which are all positive. Consequently, S > 0 and c 1 = e 1 − 2t 2 , c 2 = e 2 − 2t 1 .
If t 1 , t 2 < 0 (the hypothesis of the theorem is satisfied) we have since all terms are positive. Hence by Lemma 3.6, the heteroclinic cycle is asymptotically stable. Now suppose that t 2 < 0 < t 1 (the case t 1 < 0 < t 2 is analogous). We have since all terms are negative. By Lemma 3.6 the heteroclinic cycle is completely unstable.
The dichotomy between asymptotic stability and complete instability appears to be nongeneric for C − 2 -cycles compared to [9,Corollary 4.8]. This is due to the fact that e 2 and c 2 are not independent of t 1 . In fact, the case t 1 = 0 coincides with the degenerate situation c 2 = e 2 . Therefore, the assumption in [9, Corollary 4.8] that b 1 b 2 − a 1 + a 2 > 0 even for small t 1 > 0 cannot be satisfied here.

3.3.
Eigenvalues and eigenvectors of the transition matrix products. In the previous section we used results from [3,9] (stated as Lemmas 3.4 and 3.6) to determine the stability of the cycle. We now relate these to the hypotheses in Propositions 2.6 and 2.7 by calculating eigenvalues and eigenvectors of the transition matrix products. This is useful for our stability analysis in the higher-dimensional system in Section 4.
For α = π 2 the transition matrix product M as defined in (23) has eigenvalues λ 1 > λ 2 given by and corresponding eigenvectors If t 1 , t 2 < 0, then both eigenvalues are real and hence condition (A) is satisfied. Moreover, by the calculations in the proof of Theorem 3.7, we have so (B) is satisfied as well. Since in this case all transition matrices have only nonnegative entries, Proposition 2.6 applies and the cycle is asymptotically stable. We note that (C) is also satisfied, because 4a 1 b 1 b 2 > 0 implies that Similarly, for the components of the other eigenvector we get This is not directly related to condition (C), but will also be used in the following subsection.
On the other hand, if t 2 < 0 < t 1 , the transition matrix M 1 has a negative entry. Again by the calculations in the proof of Theorem 3.7 we have a 1 − a 2 − b 1 b 2 > 0, and therefore Thus, (C) is violated and by Proposition 2.7(a) the cycle is completely unstable. The case t 1 < 0 < t 2 is analogous for the other transition matrix product.
Note that there are other equilibria that are not part of either cycle in the heteroclinic network. For example, on SSSS all populations are phase synchronized and its stability is governed by the (quadruple) eigenvalue λ SSSS = 4r cos(2α) − 2 cos(α). For δ = 0 we have λ SSSS = λ DDSS 3 which implies that SSSS is linearly stable if the transverse eigenvalues within the corresponding subspace of each cycle are negative; cf. Section 3.2.

4.2.
Stability of the cycles. Note that by construction, the saddle SDSS has a two-dimensional unstable manifold. Hence, neither cycle can be asymptotically stable for δ ≈ 0 and α ≈ π 2 . Since the cycles are quasisimple, we can determine their stability by looking at the corresponding transition matrices. Because of the parameter symmetry, we restrict ourselves to the cycleĈ 2 ⊂ ψ 1 ψ 2 ψ 3 S without loss of generality and just write C and σ q for the remainder of this subsection.
Within the invariant subspace ψ 1 ψ 2 ψ 3 S, we have one contracting, expanding, and transverse direction with local coordinates denoted by v, w, z as above. In addition there is another transverse direction-denoted by z ⊥ in local coordinates-which is mapped to itself under the global map. The second transverse eigenvalues (those transverse to ψ 1 ψ 2 ψ 3 S) evaluate to There are two possibilities for transverse bifurcations when δ changes. If δ > 0, there is a transverse bifurcation at t ⊥ SDSS = 0. But since t ⊥ SDSS = e SDSS the other cycle of the network then ceases to exist. If δ < 0, there is a possibility of two simultaneous transverse bifurcations when t ⊥ DDSS = t ⊥ SDDS = 0. Write b ⊥ q = −t ⊥ q /e q . Again, the global maps are permutations of the local coordinate axes and the return map evaluates to In logarithmic coordinates (η, ξ, ξ ⊥ ) this gives the transition matrix that governs the stability of the cycle. Note that the upper left 2 × 2 submatrix is the same as the transition matrix (20). In order to simplify notation we write ξ 1 = SDSS and ξ 2 , . . . , ξ 6 for the subsequent equilibria of C. Assuming that we are in a parameter region where the network exists, see Theorem 4.1, we can now make the following statement about the stability of its subcycles.
Theorem 4.2. Assume that the cycle C is asymptotically stable 3 within the three-dimensional subspace it is contained in and |δ| sufficiently small. Then we have the following dichotomy.
Proof. Since t ⊥ 1 > 0 is the only positive transverse eigenvalue of an equilibrium in C, the transition matrix M 1 is the only one with a negative entry, b ⊥ 1 < 0. By Proposition 2.7 the stability of C depends on whether or not all M (q) satisfy conditions (A)-(C) in Section 2. Statement (ii) follows immediately by Proposition 2.7(a).
For (i), we want to apply Proposition 2.7(b). By [4,Lemma 3.6] it suffices to show that M (2) satisfies conditions (A) and (B), because then all M (q) satisfy (A)-(C). We calculate where µ, ν,μ,ν > 0. For a moment, suppose that δ = 0. Due to the symmetry of the system in the subspace ψ 1 ψ 2 ψ 3 S, the upper left 2 × 2 submatrix is the third power of the matrix M in (23) and we can use our calculations from Section 3.3. Note that M (2) has an eigenvalue λ = 1 with eigenvector (0, 0, 1). Its other two eigenvalues are the third powers of those of M, call them λ 1 > λ 2 , by a slight abuse of notation. Then λ max = λ 1 > 1 under the assumptions of this theorem, so conditions (A) and (B) are satisfied. Proposition 2.7(b) applies and C is f.a.s.. Since eigenvectors and eigenvalues vary continuously in δ, the same is true for |δ| sufficiently small. In order to derive expressions for the stability indices we have to find the arguments β (l) ∈ R 3 of the function F ind from Proposition 2.7. As is shown in [4], this becomes simpler if for all q = 1, . . . , 6 we have x 3 < 0 } and the convergence is demanded in every component. Clearly, this asymptotic behavior is controlled by the eigenvectors of M (q) . Consider first the case q = 2. Under our assumptions, all components of the eigenvector corresponding to the largest eigenvalue λ max > 1 have the same sign. Another eigenvector is (0, 0, 1). From Section 3.3 we know that the first two components of the remaining eigenvector have opposite signs. It follows that any x ∈ R 3 − written in the eigenbasis of M (2) must have a nonzero coefficient for the largest eigenvector. Therefore, x ∈ U −∞ (M (2) ), so (32) holds. For q = 2 note that all M (q) are similar, hence they have the same eigenvalues. Their eigenvectors are obtained by multiplying those of M (2) by M 2 , M 3 M 2 , . . . , M 6 M 5 M 4 M 3 M 2 , respectively. This involves only matrices with nonnegative entries and thus does not affect our conclusions using the signs of the entries of the eigenvectors. Therefore, (32) holds for all q = 1, . . . , 6.
Since (32) is satisfied, the only arguments β (l) ∈ R 3 that must be considered for F ind in the calculation of σ q are the rows of the (products of) transition matrices M q , M q+1 M q , M q+2 M q+1 M q and so on. Among these, we only need to take rows into account where at least one entry is negative; if there are none, the respective index is equal to +∞. Negative entries can only occur when M 1 is involved in the product, and then only in the last row. So for σ q the last row of M 1 M 6 · · · M q must be considered. Since M 2 has no negative entries and its third column is (0, 0, 1), the first two entries in the last row of M 2 M 1 M 6 · · · M q are greater than the respective entries of M 1 M 6 · · · M q , yielding a greater value for F ind . The same goes for M 3 M 2 M 1 M 6 · · · M q and so on. Thus, σ q is indeed obtained by plugging the last row of M 1 M 6 · · · M q into F ind . We get The lower right entry of all these matrices is 1, so for all q = 1, . . . , 6 we can write σ q = F ind (µ q , ν q , 1) > −∞. Since the last row of M 1 is (b ⊥ 1 , 0, 1), we have µ 1 = b ⊥ 1 and ν 1 = 0 as claimed. The recursive relations now follow immediately from (31).  Figure 3. Stability indices of the heteroclinic cyclesĈ 2 andČ 2 as in Figure 2 for (α, K, r) = ( π 2 , 0.2, 0.01). The symbol '+' for a stability index denotes 'positive and finite' and '−' denotes 'negative and finite'.
We conclude this section with a few remarks on these stability results. ) is an eigenvector of M associated with its largest eigenvalue, so both of its components have the same sign. To fulfill (C), we need sgn(u max 3 ) = sgn(u max 1/2 ). A straightforward calculation yields This condition is stronger than assuming (C), and as soon as it is satisfied, we have σ 2 = +∞.
By contrast, the indices σ 1 and σ 6 are always finite because F ind has at least one positive and at least one negative argument through b ⊥ 1 . This makes sense because they are indices along connections shared with the other cycle in the network, while σ 2 belongs to the trajectory that is furthest away from the common ones. For the other indices σ 3 , σ 4 , σ 5 there is not necessarily a negative argument, so they could be equal to +∞. From the recursive relations between the µ q and ν q we see that σ q = +∞ implies σ q−1 = +∞ for q ∈ {3, 4, 5}, which is also plausible in view of the architecture of the heteroclinic network.
Since σ q > −∞ for all q, we have shown that under the assumptions of Theorem 4.2 the cycle C is not only f.a.s., but indeed attracting a positive measure set along each of its connections. Straightforward constraints on µ q , ν q given through the definition of F ind determine the signs of all σ q and thus yield necessary and sufficient conditions for C to be even e.a.s. A simple example for such a necessary condition is b ⊥ 1 > −1, so that σ 1 > 0. This is the same as e SDSS > t ⊥ SDSS and in terms of the network parameters amounts to δ > 0, cf. Figure 3.
Similar conditions for the other σ q become increasingly cumbersome to write down explicitly and we gain little insight from them. Instead, we evaluated the stability indices (of both cycles) numerically. Two cases are illustrated in Figure 3. We conjecture that there is an open parameter  region where the assumptions of Theorem 4.2 are satisfied and the network is maximally stable (though not asymptotically stable) due to both cycles being e.a.s.. We comment further on this in the next subsection.

4.3.
Stability of the heteroclinic network. Even if the stability of all cycles that constitute a heteroclinic network is known, it is hard to make general conclusions about the stability of the network as a whole. For "simple" cases, like the Kirk and Silber network [18], a comprehensive study can be found in [19]. Based on the results in the previous section, one can draw several conclusions. If one cycle of N 2 is f.a.s.-conditions are given in Theorem 4.2-then the network itself is f.a.s. Moreover, if one cycle, sayĈ 2 , is e.a.s. and the heteroclinic trajectories inČ 2 that are not contained inĈ 2 have positive stability indices, then the network is e.a.s.-this is the case in Figure 3(b). The geometry of the two-dimensional manifold W u (SDSS) ⊂ SDψ 3 ψ 4 gives insight into the dynamics near the heteroclinic network N 2 . For simplicity, we focus on the case α = π 2 . By (28), the dynamics of the phase differences on SDψ 3 ψ 4 are given bẏ Note that if |δ| K < 2 |r| there is a (saddle) equilibrium ξ SDψ 3 D ∈ SDψ 3 D with ψ 3 = arccos(δK/2r) ∈ (0, π). For the same condition there is an analogous equilibrium ξ SDDψ 4 ∈ SDDψ 4 with ψ 4 = arccos(−δK/2r) ∈ (0, π). The stable manifolds of these saddle equilibria now organize the dynamics on SDψ 3 ψ 4 . Proof. We first give conditions on the parameters that ensure that there are no asymptotically stable sets in the invariant set (0, π) 2 ⊂ SDψ 3 ψ 4 . It suffices to show that there are no equilibria in (0, π) 2 . By direct calculation, one can verify that for α = π 2 , δ = 0 any equilibrium in (0, π) 2 must lie in {ψ := ψ 3 = ψ 4 } ⊂ (0, π) 2 . The dynamics of ψ are given byψ = (K − 4r) cos(ψ) + K. Hence, there are no equilibria if 0 < r < K/2 given K > 0; these are exactly the conditions for there to be an asymptotically stable heteroclinic cycle in each subspace by Lemma 3.1 and Theorem 3.5. Now W s (ξ SDψ 3 D ), W s (ξ SDDψ 4 ) are-as source-saddle connections-robust heteroclinic trajectories [SDSS → ξ SDψ 3 D ], [SDSS → ξ SDDψ 4 ]. These separate (0, π) 2 into three distinct sets of initial conditions which completes the proof.
The dynamics on (0, π) 2 ⊂ SDψ 3 ψ 4 are shown in Figure 4. The stable manifolds of ξ SDψ 3 D and ξ SDDψ 4 subdivide (0, π) 2 robustly into three wedges with nonempty interior that lie in the stable manifolds of SDDS, SDSD, and SDDD, respectively. In particular, this suggests that a significant part of trajectories passing by SDSS will approach SDDD which is not contained in either heteroclinic cycle of the network N 2 .
Remark 4.4. Let N be a heteroclinic network and let ξ p , p = 1, . . . , P , denote the equilibria of all its cycles. Abusing the ambiguity of Definition 2.1 4 , we call N complete [20] or clean [21] if W u (ξ p ) ⊂ P q=1 W s (ξ q ) for all p and almost complete if the set W u (ξ p ) ∩ P q=1 W s (ξ q ) is of full measure for all p and any Riemannian measure on W u (ξ p ); see also [2] for a detailed discussion in the context of coupled oscillator populations.
For N 2 to be almost complete for δ ≈ 0, the set W u (SDSS)∩(W s (SDDS)∪ W s (SDSD)) would have to be of full (Lebesgue) measure in SDψ 3 ψ 4 . However, Lemma 4.3 shows that there is a set of nonvanishing measure in SDψ 3 ψ 4 which lies in the stable manifold of SDDD, an equilibrium which is not in the network. Hence, N 2 cannot be almost complete (nor complete).
We further explored the dynamics near the heteroclinic network for M = 4 populations of N = 2 oscillators using numerical simulations with additive noise. Specifically, for (28) written asθ σ,k = ω + Y σ,k (θ)-see (2)-and independent Wiener processes W σ,k , we solved the stochastic differential equation (34)θ σ,k = ω + Y σ,k (θ) + ηW σ,k 4 Definition 2.1 is somewhat ambiguous since it does not specify how many heteroclinic connections belong to the heteroclinic network. If N2 only contains one (one-dimensional) heteroclinic trajectory (as suggested by the proof of Theorem 4.1) then it is clearly not almost complete since dim(W u (SDSS)) = 2. Strictly speaking, for the discussion of completeness we actually consider a network N2 that contains the equilibria of N2 and all connecting heteroclinic trajectories. Since the network is not stable, trajectories may go away from the heteroclinic network along the connection [SDSS → SDDD] either to return to a neighborhood of the network (top left) or to converge to the sink SSSS (top right).
using XPP [22]. As shown in Figure 5, the solutions show transitions either along the heteroclinic trajectory [SDSS → SDDS] or [SDSS → SDSD]. Since the heteroclinic network is not asymptotically stable, the dynamics also show excursions away from the network: there are trajectories that follow the heteroclinic connection [SDSS → SDDD] before either returning to the neighborhood of the network (Figure 5(a) left) or approaching the sink SSSS ( Figure 5(a) right). If the symmetry is broken by δ > 0 ( Figure 5(b)), trajectories appear to predominantly follow the principal direction [SDSS → SDDS] with the largest unstable eigenvalue as expected [23].

Discussion
Coupled populations of identical phase oscillators do not only support heteroclinic networks between sets with localized frequency synchrony but the stability properties of these networks can also be calculated explicitly. Rather than looking at dynamical systems with generic properties at the equilibria, we focus on a specific class of vector fields and obtain explicit expressions for the stability of a heteroclinic network-a feature of the global dynamics of the system-in terms of the coupling parameters. In particular, this does not exclude the possibility that stability properties depend nonmonotonously on the coupling parameters. The coupling parameters themselves can be related to physical parameters of interacting real-world oscillators, for example through phase reductions of neural oscillators [24].
Our results motivate a number of further questions and extensions, in particular in the context of the first part of the paper [2]. First, we here restricted ourselves to the quotient system; this is possible by considering nongeneric interactions between oscillator populations. The question what the dynamics look like if the resulting symmetries are broken, will be addressed in future research. Second, what happens for coupled populations with N > 2 oscillators? The existence conditions for cycles in [2] and the numerical results in [1] suggest existence of such a network, but stability conditions would rely on the explicit calculation of the stability indices [3]. In particular, the main tool used here, namely the results for quasi-simple cycles [4], ceases to apply since the unstable manifold of SDSS would be of dimension four and contain points with different isotropy; cf. [2].
How coupling structure shapes the global dynamics of a system of oscillators is a crucial question in many fields of application. Hence, our results may be of practical interest: in the neurosciences for example, some oscillators may fire at a higher frequency than others for some time before another neural population becomes more active [25]. The networks here mimic this effect to a certain extent: trajectories which move along the heteroclinic network correspond to sequential speeding up and slowing down of oscillator populations. At the same time, large scale synchrony is thought to relate to neural disfunction [26]. From this point of view, the (in)stability results of Section 4 appear interesting, since trajectories in numerical simulations may get "stuck" in the fully phase synchronized configuration SSSS. Hence, our results may eventually elucidate how to design networks that avoid transitions to a highly synchronized pathological state.

Appendix A. Phase Oscillator Populations with Nonpairwise Coupling
Interaction of phase oscillators through state-dependent phase shift may be approximated by nonpairwise coupling as shown in [1]; here we generalize these calculations to allow for arbitrary coupling topologies.
Consider a system of M populations of N phase oscillators where θ σ,k denotes the phase of oscillator k of population σ. Recall that the Kuramoto order parameter where the interaction is mediated through the coupling function (36) g(ϑ) = sin(ϑ + α) − r sin(a(ϑ + α)), a ∈ N, and the phase-shifts are linear combination of the (square of the) Kuramoto order parameters.
as an approximation for (35). Note that the interaction between different populations is through nonpairwise coupling terms: the arguments of the trigonometric functions inG (4) depend on four phase variables rather than just differences between pairs of phases.