2-Cluster Fixed-Point Analysis of Mean-Coupled Stuart-Landau Oscillators in the Center Manifold

We reduce the dynamics of an ensemble of mean-coupled Stuart-Landau oscillators close to the synchronized solution. In particular, we map the system onto the center manifold of the Benjamin-Feir instability, the bifurcation destabilizing the synchronized oscillation. Using symmetry arguments, we describe the structure of the dynamics on this center manifold up to cubic order, and derive expressions for its parameters. This allows us to investigate phenomena described by the Stuart-Landau ensemble, such as clustering and cluster singularities, in the lower-dimensional center manifold, providing further insights into the symmetry-broken dynamics of coupled oscillators. We show that cluster singularities in the Stuart-Landau ensemble correspond to vanishing quadratic terms in the center manifold dynamics. In addition, they act as organizing centers for the saddle-node bifurcations creating unbalanced cluster states as well for the transverse bifurcations altering the cluster stability. Furthermore, we show that bistability of different solutions with the same cluster-size distribution can only occur when either cluster contains at least $1/3$ of the oscillators, independent of the system parameters.

The intrinsic dimensionality of each oscillatory unit may range from d = 2 for FitzHugh-Nagumo [23] and Van der Pol oscillators [24], via d = 3 for the Oregonator [25] to d = 4 for the original Hodgkin-Huxley model [26], and even higher for more detailed physical models [27]. A system composed of N of these oscillators thus lives in a d · N -dimensional phase space, making its full investigation unfeasible even for small d and N . One can, however, circumvent this problem of increasingly large dimensions by restricting the dynamics to the center manifold of certain bifurcations. In particular, it is known that the center space of the Benjamin-Feir instability is N − 1 dimensional [20,28], and thus a reduction to the center manifold at this bifurcation allows for reducing the dimension of the problem to N − 1 and thus by a factor of ≈ d. As we show below, such a reduction lets us reveal invariant sets and bifurcation curves analytically -a difficult task in the original d · N -dimensional space.
In this work, we focus on a particular example of a globally coupled system, in which the network is composed of oscillating units called Stuart-Landau oscillators, each represented by a complex variable W k ∈ C. As opposed to phase oscillators, each Stuart-Landau oscillator has two degrees of freedom, i. e. an amplitude and a phase. With a linear global coupling, the dynamics are then given bẏ with the complex coupling constant β r + iβ i and the real parameter γ, also called the shear [29]. · indicates the ensemble mean andẆ = dW/dt. Bold face W indicates a vector containing the ensemble values [W 1 , W 2 , . . . , W N ]. For β r + iβ i = 0 the ensemble is decoupled, and each Stuart-Landau oscillator oscillates with unit amplitude and angular velocity −γ. For β r + iβ i = 0, however, a plethora of different dynamical states can be observed. These states include fully synchronized oscillations, in which all oscillators maintain an amplitude equal to one and have a mutual phase difference of zero [30], cluster states, in which the ensemble splits up into two or more sets of synchrony [31,32,14], and a variety of quasi-periodic and chaotic dynamics [33,12].
2-cluster states can be born and destroyed at saddle-node bifurcations if the number of oscillators in each cluster are different, that is, when they are unbalanced [13]. Balanced solutions with N 1 = N 2 emerge from the synchronized solution at the Benjamin-Feir instability. For N = 16 oscillators and for γ = 2, the saddle-node bifurcations for different unbalanced cluster distributions N 1 = N 2 and the Benjamin-Feir instability are depicted in Fig. 1, as a function of the coupling parameters β r and β i . Here, all the 2-cluster solutions exist locally in parameter space below their respective saddle-node bifurcation curve, that is for smaller β r values. Up to the Benjamin-Feir instability they coexist with the stable synchronized solution. Descending from large β r values, notice that the most-unbalanced cluster state with N 1 : N 2 = 1 : 15 is created first. The more balanced cluster states are born subsequently, depending on their distribution, until eventually the balanced cluster state N 1 : N 2 = 8 : 8 is born at the Benjamin-Feir instability. At β r = −(1 − √ 3γ)/2, β i = (−γ − √ 3)/2, there exists a codimension-two point where the saddle-node bifurcations of all cluster distributions coincide. This point is called a cluster singularity [32]. Note that the qualitative picture in Fig. 1 does not change when increasing the total number of oscillators N . For large numbers N → ∞ of infinitely many oscillators we expect a bow-tie-shaped band of saddle-node bifurcation curves, ranging from the saddle-node bifurcation of the most unbalanced cluster state to the Benjamin-Feir instability. As argued in Ref. [32], the cluster 2-Cluster Fixed-Point Analysis of Mean-Coupled Stuart-Landau Oscillators in the Center Manifold singularity can thus be viewed as an organizing center. By projecting the dynamics close to the Benjamin-Feir instability onto its center manifold, we aim to obtain further insights into the properties of this organizing center, and to elucidate the clustering behavior near it.  Figure 1: The Benjamin-Feir instability involving the 8 : 8 cluster (dark blue) and the different saddle-node curves creating the unbalanced cluster solutions, N 1 = N 2 , in the β i , β r plane with γ = 2 and N = 16. Each curve belongs to a particular cluster distribution N 1 : N 2 , and is obtained with numerical continuation using auto-07p [34,35]. Note the position of the cluster singularity at The remainder of this article is organized as follows: In Sec. 2, we pass to a corotating frame and introduce the average amplitude R, the deviations from the average amplitude r k , the deviations from the mean phase ϕ k . Using this corotating system, we discuss how one can describe the dynamics in the center manifold, see Sec. 3. In Sec. 4 we derive the parameters for the dynamics of x k . Detailed calculations are provided in Appendix C, for convenience. Based on the parameters in the center manifold, we study the bifurcations of 2-cluster states and the role of the cluster singularity in the center manifold, in Sec. 5. We conclude with a detailed discussion of our results and an outlook of future work. For a detailed mathematical analysis of the dynamics of 2-cluster states in the center manifold, see the companion paper Ref. [36].

Variable transformation into corotating frame
Notice that Eq. (2) is invariant under a rotation in the complex plane W k → W k exp iφ. This invariance can be eliminated by choosing variables in a corotating frame, thus effectively reducing the dimensions of the system from 2N to 2N − 1.

2-Cluster Fixed-Point Analysis of Mean-Coupled Stuart-Landau Oscillators in the Center Manifold
In particular, we express the complex variables W k in log-polar coordinates W k = exp(R k + iΦ k ). Then Eq. (2) turns intoṘ with k = 1, . . . , N −1, the abbreviations shown in Tab. 1 and the new coordinates summarized in Tab. 2 (see Appendix A for a derivation). Hereby, · symbolizes the deviation from the ensemble mean · , and R and Φ are the ensemble mean Table 1: Abbreviations logarithmic amplitude and phase, respectively. The logarithmic amplitude and phase deviation of each oscillator from their averages are r k and ϕ k . Notice that through this construction, the averages of these deviations vanish. Furthermore, bold face of a variable, e.g. x, symbolizes the set of the respective ensemble variables {x 1 , x 2 , . . . , x N }. To simplify notation, r k + iϕ k is abbreviated by the complex variable z k . The transformation into Eqs. (3a) to (3c) has the advantage that the resulting equations are independent of the mean phase Φ. A change of Φ corresponds to a uniform phase shift of the whole ensemble in the complex plane, which in turn means that periodic orbits in the Stuart-Landau ensemble, Eq. (2), correspond to stationary solutions in the transformed system, Eqs. (3a) to (3c). Thus, we can ignore the mean phase Φ in our subsequent analysis. Synchronized oscillations correspond to R = 0, r k = 0, Φ = −γt and ϕ k = 0. The stability of this equilibrium can be investigated using the eigenspectrum of the Jacobian evaluated at this point. Due to the S N -symmetry of the solution and the S N -equivariance of the governing equations, the Jacobian becomes block-diagonal, and thus has a degenerate eigenvalue spectrum [21,15], see Appendix B: • There is one singleton eigenvalue λ 1 = −2 < 0, corresponding to an eigendirection affecting all oscillators identically. That is, this direction v 1 shifts the amplitude of the synchronized motion but does not alter its symmetry.
• There is the eigenvalue λ + = −1 − β r + 1 − β 2 i − 2β i γ =: −1 − β r + d which becomes zero at the Benjamin-Feir instability and is of geometric multiplicity N − 1. The corresponding directions correspond to 2-cluster states, with each direction corresponding to one cluster distribution N 1 : N 2 . Up to conjugacy, we arrange here the units such that the first N 1 oscillators correspond to the same cluster. All 2-clusters with the same distribution but different assignments of the oscillators then belong to the same conjugacy class.
• Finally, there is the eigenvalue λ − = −1 − β r − d which is negative close to the synchronized solution, which has a geometric multiplicity of N − 1 and whose eigendirections also have S N1 × S N2 -symmetry.

Center manifold reduction
In the following, we calculate an expansion to third order of the (N −1)-dimensional center manifold which corresponds to the Benjamin-Feir instability at λ + = 0 = −1 − β r + d. In order to do so, it is useful to introduce the coordinates such that Here we use the notations γ = 2γ + β i and d as defined above. See Appendix B for a derivation. The variables x k describe the dynamics in the (N − 1)-dimensional center manifold tangent to y k = 0 ∀k, while y k together with R describe the dynamics in the stable manifold tangent to x k = 0 ∀k. Note that the center-manifold must be S N -invariant. In addition, the global restrictions r = ϕ = 0 and thus x = y = 0 must hold. Therefore, the general form of the center manifold up to quadratic order must follow with the coefficients a = a (β i , γ) and b = b (β i , γ). Here, we use tangency of our coordinates R and y k , that is, Since the Benjamin-Feir instability β r = d − 1 is of codimension one, the three-dimensional parameter space (β r , β i , γ) becomes two-dimensional. The parameters in the center manifold thus only depend on β i and γ. By S N -equivariance, the reduced dynamicsẋ k in the center manifold, up to cubic order, must be of the formẋ see also Refs. [10,20], with the parameters A = A (β i , γ) and B = B (β i , γ) and C = C (β i , γ).

Derivation of the parameters a, b, A, B and C
In this section, we discuss the approach to calculate the coefficients a, b, A, B and C for the dynamics in the center manifold. See Appendix C for complete details. First, we determine b. In particular we observe thaṫ holds. Since λ + = 0 at the bifurcation,Ṙ up to second order in x k must vanish. Therefore, expressing z k = r k + iϕ k and r k , ϕ k in terms of x k in Eq. (3a), we can compute b by comparing the coefficients of the x 2 : the terms in front of x 2 must thereby vanish. This allows us to with γ and d as defined above. Analogously, we can calculate a using Eqs. 3a and 3b up to second order in x k and employinġ This means we can use 2dẏ k =ṙ k + (d − 1)/γ φ k , substitute the z k with x k in Eqs. 3a and 3b and keep terms up to O x 2 k . Comparing the coefficients in front of x 2 k then results in Finally, we can calculate A, B and C using Taking Eqs. 3b and 3c and the coefficients a and b obtained above, we can evaluate this equality up to cubic order, yielding the coefficients Together with λ + , the expressions for A, B and C fully specify the dynamics in the center manifold based on the original parameters γ, β r and β i . By rescaling time and x k in Eq. (10), the number of independent parameters can be reduced to two, see Ref. [36]. For simplicity, we use the unscaled equation as in Eq. (10) here.

Clustering and cluster singularities in the center manifold
As shown in Fig. 1 for N = 16 oscillators, we observe a range of saddle-node bifurcations creating the different 2-cluster states. The expressions for λ + , A, B and C above determine the corresponding parameter values in the center manifold. The respective λ + and A values for the numerical curves shown in Fig. 1 are depicted in Fig. 2 as dashed curves. Notice that the Benjamin-Feir curve corresponds to the line λ + = 0. Furthermore, we can derive the saddle-node curves creating unbalanced 2-cluster states in the center manifold analytically, see Appendix D. In particular, for unbalanced cluster solutions, with α = N 1 /N 2 . The respective analytical curves for N = 16 are shown as solid curves in Fig. 2. Notice the close correspondence between the mapped bifurcation curves from the full system and the bifurcation curves determined in the center manifold. For more balanced solutions, the saddle-node curves obtained from the Stuart-Landau ensemble depart more strongly from the saddle-node curves calculated analytically in the center manifold. We expect this to be due to the cubic truncation of the flow in the center manifold, thus limiting its accuracy away from the Benjamin-Feir curve.
Note that to obtain the curves in Fig. 2, we fix γ = 2 and vary β i , β r . We then use the expressions for A(β i , β r ), B(β i , β r ) and C(β i , β r ) to get the parameters in the center manifold. Thus the parameters A, B and C lie on a two-dimensional manifold. For the curves shown in Fig. 2, we furthermore use Eq. (16), yielding one-dimensional curves. The curves are, however, not exactly parabolas, since B and C vary in addition to A, which is not shown in Fig. 2. For all subsequent figures, we use the values of C = −1 and B = −2/(2 √ 3 − 3) at the cluster singularity for γ = 2, which can be obtained analytically. See Ref. [36] p. 36 for a derivation.
Furthermore, from Fig. 2 we observe that A = 0, in addition to λ + = 0, at the cluster singularity. This means that this codimension-two point is distinguished by vanishing quadratic dynamics in the center manifold, cf. Eq. 10. In addition, it serves as an organizing center for the saddle-node bifurcations of the unbalanced cluster states: At the saddle-node bifurcation, we have in the center manifold for a cluster state with B < 0 and C < 0 for the range of β i , β r considered here (not shown), see Appendix D. This means that for negative A values, the saddle-node curves occur at positive x 1 , for positive A values at negative x 1 , and for A = 0, at the cluster singularity, all saddle-node bifurcations occur at the synchronized solution x k = 0. This behavior can indeed be observed in the Stuart-Landau ensemble, see Fig. 6 of Ref. [32].  Fig. 1 using the expressions for λ + (β r , β i , γ) = −1 − β r + d, and A(β r , β i , γ), cf. Eq. (13). The solid curves λ sn indicate the saddle-node bifurcations of the unbalanced cluster states obtained analytically in the center manifold, see Eq. (16). Note that close to the cluster singularity, where analytical expansions work best, numerical continuation fails due to the concentration of solutions in phase space.
The unbalanced cluster states do, in general, not emerge as stable states from the saddle-node bifurcations. Rather, one of the two branches created at the saddle-node bifurcation is subsequently stabilized through transverse bifurcations involving 3-cluster solutions with symmetry S N1 × S N2 × S N3 , also called secondary branches [10]. For a more detailed discussion on secondary branches, see also Refs. [37,28] In order to explain this in more detail, we follow Ref. [37] Section 4. Note that each N 1 : N 2 2-cluster solution is invariant under the action of the group S N1 × S N2 . From this, it follows that one can block-diagonalise the Jacobian at the 2-cluster solutions S N1 × S N2 . In doing so, one can calculate the (N 1 − 1)-degenerate eigenvalue µ 1 describing the intrinsic stability of cluster Ξ 1 , that is its stability against transverse perturbations. Note, however, that a cluster of size 1 cannot be broken up. Following Ref. [10] p. 23 and using isotopic decomposition, the eigenvalue µ 1 can be expressed as Here, J ij | Ξ1 denotes ∂f i /∂x j , with the respective x i and x j in cluster Ξ 1 and f i being the right hand side of Eq. (10). Without loss of generality, we assume in the following that Ξ 1 is the cluster with the smaller number of oscillators, that is, N 1 ≤ N 2 or α ≤ 1. Evaluating the Jacobian, one obtains that the eigenvalue µ 1 changes sign at 2-Cluster Fixed-Point Analysis of Mean-Coupled Stuart-Landau Oscillators in the Center Manifold Analogously, the transverse stability of cluster Ξ 2 is described by which changes sign at Hereby, µ 2 describes the intrinsic stability of cluster Ξ 2 . Furthermore notice that for the balanced cluster α = 1 and therefore λ +,1 = λ +,2 . Since both clusters contain an equal number of units, their respective intrinsic stabilities change simultaneously. The saddle-node curve creating the 4 : 12 cluster is shown as solid orange curve. The Benjamin-Feir line is shown in blue, with the λ +,1 = λ +,2 curve for the balanced 8 : 8 cluster state depicted as dotted blue curve. The 4 : 12 cluster is stable in the two regions between the respective λ +,1 = 0 and λ +,2 = 0 curve. The balanced cluster state is stable above the dotted blue curve.
In Fig. 3, λ sn , λ +,1 and λ +,2 are shown as solid, dotted and dash-dotted orange curves, respectively, for the 4 : 12 2-cluster state. The Benjamin-Feir instability, where the balanced cluster state is born, is drawn as a solid blue line at λ + = 0, and the transverse bifurcation curve λ +,1 = λ +,2 , where the balanced cluster state is stabilized, is drawn as a dotted blue curve. See Fig. 4 for the respective curves for a range of cluster distributions. Fig. 3 can be interpreted as follows: Coming from negative λ + values, the unbalanced 4 : 12 cluster state is born at λ sn (4 : 12) (solid orange). However, this 2-cluster state is unstable for the parameter values considered here: the cluster Ξ 1 with 4 units is intrinsically unstable with µ 1 > 0 and µ 2 < 0. At the dotted orange curve, µ 1 changes sign, rendering the 4 : 12 cluster state stable. Subsequently, at the dash-dotted orange curve, µ 2 changes sign, leaving the cluster Ξ 2 with 12 units intrinsically unstable and thus the 4 : 12 cluster unstable.

2-Cluster Fixed-Point Analysis of Mean-Coupled Stuart-Landau Oscillators in the Center Manifold
The qualitatively same behavior can be observed for any cluster distribution, cf. Fig. 4, except for the most unbalanced state (1 : 15). There, cluster Ξ 1 cannot be intrinsically unstable, since it contains only one unit. This means that this cluster solution is born stable in its saddle-node bifurcation, and becomes unstable only at λ 2 when µ 2 = 0. See also the bottom right plot in Fig. 4. In particular, λ sn = A 2 /4B for α = 0, see Eq. (16), coincides with λ +,1 = A 2 /4B for α = 0, cf. Eq. (17). Furthermore, it is worth noting that the stable patches in parameter space overlap for different cluster distributions. This means that there is a multistability of different stable 2-cluster states. Notice that these results are in close correspondence with the behavior observed in the full Stuart-Landau ensemble, compare, for example, Fig. 3 with Figs. 4b and 5b in Ref. [32]. λ sn and λ +,1 are continuous functions of α. For N → ∞, this means that there are continuous bands of bifurcation curves: Going from λ sn (α = 0) = A 2 /4B to λ sn (α = 1) = 0, there is band of saddle-node bifurcations creating the unbalanced cluster solutions. This band becomes infinitesimally thin at the cluster singularity A = 0, giving it a bow-tie like shape. From λ +,1 (α = 0) = A 2 /4B to λ +,1 (α = 1) = (−B − C)A 2 /B 2 , the transverse bifurcations of the smaller cluster stretch from the saddle-node curve of the most unbalanced cluster state to the transverse bifurcations of the balanced cluster state where λ +,1 is maximal, again yielding a bow-tie like shape in the A, λ + plane. Since λ +,2 has a pole at α = 1/2, the interpretation is a bit more involved. First, for the balanced cluster state α = 1: and thus λ +,2 and λ +,1 coincide. For the most unbalanced solution α = 0: λ +,2 (α = 0) = 0. This means the larger cluster of the most unbalanced solution becomes unstable exactly when the balanced solution is born, that is, at the Benjamin-Feir instability λ + = 0. For intermediate α values, however, the λ +,2 curve becomes steeper and infinitely steep at α = 1/2, with the tip reaching to the cluster singularity. This can also be observed in Fig. 4, where the parabola becomes thinner when going from the 6 : 10 to the 5 : 11 cluster states, and subsequently broadens again until the 1 : 15 cluster. Altogether, the λ +,2 curves fill out the half plane λ + ≥ 0 except the line A = 0. These three bow-tie like regions of λ sn , λ +,1 and λ +,2 become infinitesimally thin and thus singular only at the cluster singularity λ + = 0, A = 0.
The bifurcation scenario can be better visualized by plotting the λ sn , λ +,1 and λ +,2 as a function of the cluster size N 1 /N , see Fig. 5. It depicts the λ + values of the saddle-node bifurcations creating the 2-cluster states (λ sn , blue) and the two transverse bifurcations (Eqs. (17) and (18)) altering the stability of the 2-clusters, with λ +,1 in green and λ +,2 in orange.
When increasing λ + coming from negative values, all cluster states with N 1 /N = 1/2 are born in the saddle-node bifurcation λ sn . Note that in fact two solutions for each N 1 /N are created this way. In Fig. 5, one can observe that for the most unbalanced state N 1 /N → 0, the transverse bifurcation stabilizing the smaller cluster λ +,1 occurs immediately after the saddle-node bifurcation creating that cluster. This bifurcation alters the stability of one of the two solutions born in the saddle-node bifurcation, and in particular renders the smaller of the two clusters in that solution stable to transverse perturbations. For the parameter regime considered here (A = −0.2, B = −2/(2 √ 3 − 3) and C = −1), this solution is in fact stabilized at this bifurcation, that is for λ + > λ +,1 . For N 1 /N < 1/3, the respective 2-cluster solution remains stable until λ +,2 , where the larger cluster becomes unstable, thus rendering the whole solution unstable. This can, for example, be observed for the 4 : 12 cluster-size distribution, see Fig. 6(top). There, the variable of one cluster, x 1 , is plotted as a function of the bifurcation parameter λ + . The blue dot on the left marks the saddle-node bifurcation wherein the two 4 : 12 solutions are created. Initially, both solutions are unstable. At λ +,1 (orange dot), one of them is stabilized and at λ +,2 (green dot), it is subsequently destabilized.
For N 1 /N > 1/3, the scenario is different. There the solution that got stabilized at λ +,1 remains stable for all λ + > λ +,1 . The bifurcation λ +,2 instead occurs at the second cluster solution created at the saddle-node bifurcation. This is illustrated more clearly in Fig. 6(bottom) for the 7 : 9 cluster solution. One of the two solutions becomes stable at λ +,1 , marked by an orange dot and as discussed above. Since N 1 /N = 7/16 > 1/3, this solution remains stable for all λ + > λ +,1 . The second solution (upper curve in the bottom part of Fig. 6) first passes the synchronized solution at the Benjamin-Feir bifurcation λ + = 0 and finally becomes stabilized at λ +,2 , marked by a green dot. λ +,2 diverges at the pole N 1 /N = 1/3, separating the two scenarios shown in Fig. 6. There the bifurcation switches from the solution with negative x 1 (which, for λ + → ∞, diverges to −∞) to the solution with positive x 1 (which, for λ + → ∞, diverges to +∞).
Notice how for the cluster distribution N 1 /N = 7/16 the two 2-cluster solutions are bistable for λ + > λ +,2 . That is, there exist two stable 2-cluster solutions with different x 1 but the same cluster size ratio 7 : 9 that are both stable. This, in fact, has also been observed in the Stuart-Landau ensemble, see for example Fig. 6 in Ref. [32]. Note that the singularity of λ +,2 at N 1 /N = 1/3 (α = 1/2) is independent of the parameters A, B and C, see Eq. (18). This means that bistable solutions created as described above can in general only exist for N 1 /N > 1/3. The saddle-node curves creating the unbalanced cluster solutions are represented as solid curves, which correspond to the shaded curves in Fig. 2 with the same color coding. The Benjamin-Feir line is shown in blue. The unbalanced cluster states are stable above the respective dotted curve and below the dash-dotted curve, except for the 1 : 15 cluster, which is stable already at the saddle-node bifurcation. For the 2 : 14 cluster, the dotted and solid curves do not coincide but lie very close in parameter space.

Conclusion and Outlook
In this paper, we showed how one can map a system of globally coupled Stuart-Landau oscillators onto the (N − 1)dimensional center manifold at the Benjamin-Feir instability. Thereby, we observed that the bifurcation curves at which 2-cluster solutions are born closely resemble their counterparts in the original oscillatory system. This allowed us to investigate a codimension-two point called cluster singularity, from which all these bifurcation curves emanate. In the center manifold, we saw that this point corresponds to a vanishing coefficient A = 0 in front of the quadratic term of the equations of motion. Due to the reduced dynamics in this manifold, we were able to obtain stability boundaries for 2-cluster states analytically. This allows for the more detailed investigation of the bow-tie-shaped cascade of transverse bifurcations that govern the stability of these 2-cluster states, highlighting the role of the cluster singularity as an organizing center. The observed behavior is hereby independent of the oscillatory nature of each Stuart-Landau oscillator, but a result of the S N -equivariance of the full system. These findings may thus facilitate our understanding of this codimension-two point, and of clustering in general, even beyond oscillatory ensembles. Through this reduction to the center manifold, we could calculate the bifurcation curves creating the cluster solutions (λ sn ) and altering their stability (λ +,1 and λ +,2 ) analytically. This allowed us to investigate when stable 2-cluster solutions exist more systematically, and in particular revealed when different solutions with the same cluster-size distribution are bistable (cf. Fig. 6). The relative cluster size N 1 /N = 1/3 seems to be a general lower limit for such a bistable behavior. The bifurcation scenario of how states with different cluster size ratios N 1 /N are created is thereby different from the Eckhaus instability [38] in reaction-diffusion systems. There, solutions of different wavelengths are created through supercritial pitchfork bifurcations at the trivial solution and subsequently stabilized through a sequence of subcritical pitchfork bifurcations involving mixed-mode states. In our case, the different 2-cluster states are created in saddle-node bifurcations and stabilized at λ +,1 at a single equivariant bifurcation point involving 3-cluster states. However, the detailed interaction between 2-and 3-cluster states still remains an open topic for future research. Note that the cubic truncation of the flow in the center manifold has a gradient structure [36]. This means that we can assign an abstract potential to each of the cluster distributions for a particular set of parameters λ + , A, B and C. Is there a particular cluster distribution with a minimal potential value? What is its role in the dynamics between these cluster distributions? The companion paper [36] addresses some of these dynamical questions. Here, we fixed the parameter γ = 2 in the full Stuart-Landau system, and varied the coupling parameters β r , β i . This restricts our analysis to a small region in parameter space. It is important to mention that for different parameter regimes, a qualitatively different behavior close to the cluster singularity might be observed [36]. As discussed in Sec. 2, the Stuart-Landau ensemble permits the transformation into a corotating frame. This turns limit-cycle dynamics into fixed-point dynamics and thus greatly facilitates the reduction onto the center manifold. For more general oscillatory ensembles, such as systems composed of van der Pol or Hogdkin-Huxley type units, the transformation to a corotating frame may be more cumbersome or not even possible. If the coupling between such units is of a global nature, we expect, however, that the nesting of bifurcation curves creating different cluster distributions, cf. Fig. 2, can also be observed in these systems.

Acknowledgement
FPK thanks BF for the hospitality and the exciting discussions at the Freie Universität Berlin. BF gratefully acknowledges the deep inspiration by, and hospitality of, his coauthors at München who initiated this work. This work has also been supported by the Deutsche Forschungsgemeinschaft, SFB910, project A4 "Spatio-Temporal Patterns: Control, Delays, and Design", and by KR1189/18 "Chimera States and Beyond".

B Linearization
Linearizing the dynamics of the transformed system, Eqs. (3a) to (3c), at the equilibrium R = 0, r k = ϕ k = 0, z k = 0, and using the fact that r = 0, z = 0, see Tab. 2, one gets The Jacobian thus has the eigenvalues • Eigenvalue λ 1 = −2 with eigenvector v 1 = 1, 0, 0 . and two eigenvalues of geometric multiplicity N − 1 given by the eigendecomposition which gives Here, we assume 1 − β 2 i − 2β i γ > 0, that is real λ ± . For an analysis of the case 1 − β 2 i − 2β i γ < 0, see Ref. [39]. The eigenvectors corresponding to these two eigenvalues can be obtained using  For λ + , one thus obtains Choosing thus solving the equality above. The constraint Eq. (A.1), together with r = ϕ = 0, defines an Choosing 2) solves the conditions above. In particular, The constraint Eq. (A.2), together with r = ϕ = 0 define an (N − 1)-dimensional subspace of R 2N −1 . Now, one can define the eigencoordinates x k describing the dynamics in the space defined by the constraint Eq. (A.1), the center space of the bifurcation, and eigencoordinates y k , describing the dynamics in the space defined by the constraint Eq. (A.2). These two sets of variables, together with R, can then be used to describe the full system.

C Parameter Derivation
In this section of the appendix, we derive expressions for the parameters a, b, A, B and C as a function of the parameters γ, β r and β i from the Stuart-Landau ensemble. Hereby, we will use the condition that R and the y k are tangential, that is, d dx k R x=0 = 0 and d dx k y k x=0 = 0.