Design principle of multi-cluster and desynchronized states in oscillatory media via nonlinear global feedback

A theoretical framework is developed for the precise control of spatial patterns in oscillatory media using nonlinear global feedback, where a proper form of the feedback function corresponding to a specific pattern is predicted through the analysis of a phase diffusion equation with global coupling. In particular, feedback functions that generate the following spatial patterns are analytically given: (i) 2-cluster states with an arbitrary population ratio, (ii) equally populated multi-cluster states and (iii) a desynchronized state. Our method is demonstrated numerically by using the Brusselator model in the oscillatory regime. Experimental realization is also discussed.


Introduction
Feedback control is a powerful method of regulating spatio-temporal dynamics and has been studied in a wide variety of fields including physics, chemistry, biology and medical science [1]. For example, formation of various clustering patterns has been realized in the Belousov-Zhabotinsky reaction [2,3] and in the catalytic CO oxidation reaction on Pt [4]- [6]. The catalytic CO oxidation systems have also been studied for the suppression of chemical turbulence [5,7]. Moreover, considerable attention has been paid to feedback devices that suppress the pathological synchronization in the brains of Parkinson's disease patients [8]- [13].
In many cases systems to be controlled are spatially extended, and reaction-diffusion systems provide a good model for the study of pattern controlling. Theoretical analyses based on reaction-diffusion systems have been done for the Belousov-Zhabotinsky reaction [3,14,15] and CO oxidation [16]- [18]. However, so far, only empirical control has been achieved for such spatially extended systems, including the above-mentioned pioneering experimental works [2,3,6]; there has been no general theory that quantitatively relates feedback inputs to spatial patterns.
On the other hand, for discrete oscillator systems, such a quantitative feedback control methodology has been established very recently by Kiss et al [19] and Kori et al [20]. Their method is based on a phase model described by 3 model, and any coupling function can be designed by applying an appropriately constructed feedback signal to a population of oscillators. Hence the population of oscillators can be steered to a desired synchronization behavior by taking the following two steps: (i) find a coupling function that results in a desired synchronization behavior in (1) and (ii) construct an appropriate feedback signal that yields the coupling function. A major advantage of their methodology is that the phase model can be constructed from experimentally measurable quantities only; detailed information on the intrinsic dynamics of the system is not necessary. The validity and robustness of their methodology have been confirmed both experimentally by using electro-chemical oscillators [19,20] and numerically [20]. In this paper, by utilizing the above methodology of Kiss, Kori, Hudson and Rusin, we develop a general theory for the global feedback control of spatially extended oscillatory media. Our approach is also based on a phase model. Since the existence of diffusive coupling plays a crucial role in the development of spatial patterns in oscillatory media, our phase model inevitably includes both diffusive and global coupling, in contrast to discrete oscillators. By studying such a phase model, we find coupling functions leading to the following spatial patterns characterized by the distribution of phases: (i) 2-cluster states with specified population ratios, (ii) equally populated multi-cluster states and (iii) a desynchronized state. Moreover, we propose a new nonlinear feedback function without time delay, which is more convenient for designing various coupling functions than that used in previous work [19,20]. We numerically demonstrate our proposed method by using a particular reaction-diffusion model and reproduce all three patterns mentioned above with theoretically predicted feedback parameters. This paper is organized as follows. In section 2, we present the basic ideas of our control methodology for oscillatory media in detail. In section 3, we give a detailed analysis of the phase diffusion equation with special coupling functions that yield the three spatial patterns mentioned above. The numerical demonstration of the theory by using the Brusselator is given in section 4. Experimental realization is discussed in section 5.

General control methodology
Our approach for the control of oscillatory media is closely related to the method recently proposed for the population of oscillators [19,20]. The dynamics of discrete, identical limitcycle oscillators under global feedback are described by the following nonlinear dynamical equations: where u i is the state vector of the ith oscillator (i = 1, . . . , N ), F is a nonlinear function describing a limit-cycle oscillation, K is the coupling strength, h(u i ) represents the feedback and e is a unit vector with only one nonzero component: we have assumed that the feedback is additively applied to the system. When the coupling is weak, by treating the second term as a small perturbation the system is reduced to the phase model (1) [21]. In this phase description (1), synchronization behavior depends solely on the coupling function (φ), and therefore one can control the synchronization behavior of the system if the coupling function is freely given. It has been shown [19,20] that this can be done by applying a properly designed external feedback signal h(u i ) to the 4 oscillator system. This method relies on the fact that the coupling function is the convolution of the feedback h and the phase response function Z (φ), which characterizes the sensitivity of the phase to a weak external perturbation (see appendix A).
Since oscillatory media can be regarded as a population of oscillators that are diffusively connected, we argue that the same method works for shaping the coupling function in the phase description of oscillatory media. Consider a d-dimensional reaction-diffusion system with a global coupling: where u(x, t) is the state vector,D denotes the diffusion matrix, F(u) is a reaction term that generates a limit-cycle oscillation with the frequency ω, K is the coupling strength, e is the same as above and h(u) represents the feedback integrated over the entire space S. We assume that the system is Benjamin-Feir stable, i.e. the system undergoes spatially uniform oscillation when K = 0 (external control is absent). Following the standard procedure (see appendix A), we obtain where φ(x, t) is the phase of the local oscillation, and α and β are constants determined by the property of the oscillatory medium. As in the case of discrete oscillators, the coupling function (φ) can be arbitrarily shaped by properly designing the feedback signal h. For discrete oscillators, n-cluster states can be generated from the coupling function that contains nth harmonics [22]. Even in the case of oscillatory media, the coupling function is expected to work in the same way to stabilize the clustering pattern. A distinct problem here, however, is that the spatial patterns are not solely determined from the global coupling but from the interplay between the diffusive coupling and the global coupling, which makes the analysis much more complicated. For example, when a clustering pattern forms, interfaces appear between the clusters due to the diffusive coupling. Then, controlling the interface motion needs to be done to obtain the desired clustering patterns.
Hence we take the following strategy. We start from a coupling function with which discrete oscillators described by equation (1) exhibit clustering or desynchronization. We then study the phase diffusion equation (4) with this coupling function to find the resulting spatial patterns. Once the relation between the pattern and the coupling function is obtained, the corresponding feedback function can be found by following the same procedure as the discrete oscillators case. The key to carrying out this strategy is to find proper coupling functions that allow for analytical treatment of the phase equation. Although analytical treatment is easy for a simple coupling function as (φ) = sin(φ + θ) with a parameter θ [1], only poor spatial patterns appear with such a coupling function; higher harmonics in the coupling function are responsible for the formation of complex spatial patterns, including phase clustering behavior. In the next section, we propose such analytically tractable coupling functions that produce clustering states and the desynchronized state.

Analysis of the phase model
In this section, we study a one-dimensional phase diffusion equation with a global coupling (4). Here, we propose coupling functions that yield interesting spatial patterns. We are especially interested in 2-cluster states with specified population ratios, equally populated multi-cluster states and the desynchronized state.

Two-cluster states: numerical investigation
Here we focus on the 2-cluster states with an arbitrary population ratio. In particular, we look for well-defined 2-clusters, with the phases maximally separated by π .
Discrete oscillators are known to form various clustering states, the behavior is entirely governed by the form of the coupling function [22]. In particular, the following coupling function yields 2-cluster states with the phase difference equal to π: where γ and θ are parameters (see appendix B). Equation (1) with this coupling function has a family of 2-cluster solutions with different population ratios of the two clusters, which are stable for some range of population ratios. Thus, starting from a random initial condition, the system converges to a 2-cluster state with its population ratio determined from the initial condition.
Using this coupling function, we numerically investigate the phase diffusion equation (4) in one dimension with the system size S = L, taking θ as a control parameter and γ = 0.3. This equation is solved with the flux-free boundary condition by using the second-order Euler scheme, with the spatial and time intervals being set as x = 0.1 and t = 0.01, respectively. We set L = 100, K = 0.1 and α = 0.384 × 10 −2 . We assign several values to β to make the phase diagram given below. Otherwise we set β = 1.089 × 10 −2 . This special choice of α and β is for later comparison to the Brusselator model. Figure 1 shows the phase diagram obtained by varying θ for each δ ≡ β/α with several values of β and a fixed value of α. As expected, there exists a finite range of stationary 2-cluster states. Note that, as opposed to discrete oscillators, here the population ratio between the two clusters is uniquely determined for fixed δ and θ . Increase (decrease) in θ widens the phaseadvanced (retarded) region. At some critical value of θ the stationary state becomes unstable, leading to the recurrent 2-cluster state, where the following process occurs in a repeated way (see figure 2(a)): after a long transient of a quasi-stationary 2-cluster state, a new cluster sprouts out of the phase-advanced (retarded) cluster. Then the two interfaces propagate and one of the clusters disappears, the system returning to the 2-cluster state. Such dynamics have been reported in the CO oxidation model [16], although investigated only numerically.
To characterize the patterns, we introduce the lth order parameters (l = 1, 2, . . .): For 2-cluster states, |σ 1 | indicates an approximate population disparity between the two clusters and 1 − |σ 2 | indicates the ratio of the interface width to the system size L. Note the two different timescales in figure 2(b), each corresponding to the emergence of a new cluster and a slow drift of the interface. As θ exceeds the threshold around ±π/2, the recurrent 2-clusters turn into the desynchronized state (figure 1(c)), where σ 2 almost vanishes.

Two-cluster states: analytical investigation
Here, we analytically investigate the 2-cluster states numerically found above. The analysis can be done by taking the L → ∞ limit. We derive analytical forms of σ 1 and σ 2 as functions of θ and give the stability boundaries of the stationary 2-cluster shown in figure 1. We move to a co-rotating frame so that the phase of the phase-retarded cluster is fixed at φ = 0.
As we can see from the numerical result, the profile of a 2-cluster state can be decomposed into three regions: the phase-retarded cluster denoted by A (φ = 0), the phase-advanced cluster denoted by B, and the interface. Also, from the numerical observation, it is implied that the instability leading to the recurrent 2-clusters appears from the clustered region, while the interface remains stable. Hence in the analysis below, we assume that the interface does not contribute to the stability. This separation of the regions becomes well defined for large L. When the interface width is negligible compared to the system size, the two order parameters become real. In particular, in the steady state, we have σ 2 = 1, so that σ 1 is the only relevant order parameter.
Consider the dynamics of the cluster A. The contribution of the interface comes from the global coupling represented as the integral in equation (4). Since the interface width is vanishingly small, the interface region itself does not affect the dynamics. The remaining effect of the interface comes indirectly through the interface motion that varies the population ratio of the two clusters. However, since the timescale of the interface motion is O(1/L), as shown below, the population ratio can be treated as constant. Thus in this limit the dynamics of the clusters are independent of the interface motion. When the interface can be negligible, equation (4) has the solutions φ(x) = 0 for x ∈ A and φ(x) = π for x ∈ B, where the population ratio is arbitrarily given.  The stability analysis can be performed in the same way as for discrete oscillators (see appendix B). The only difference is the contribution from the diffusive coupling, which turns out to be negligible in the large L limit. Two modes of fluctuation occur in the 2-cluster state: intercluster and intra-cluster fluctuations. The inter-cluster mode is a fluctuation of the phase between the clusters, with each cluster oscillating uniformly. The eigenvalue associated with this mode is given by λ inter = −1 − 2γ cos θ . Thus, by choosing |γ | < 1 2 , we can keep this mode stable. On the other hand, for the intra-cluster mode, a fluctuation within a cluster can be unstable. The eigenvalues associated with the clusters A and B with the wavenumber k are given by λ (A) intra = −αk 2 + 2 p − 1 − 2γ cos θ and λ (B) intra = −αk 2 + 1 − 2 p − 2γ cos θ , respectively, where p, the population ratio, is the area fraction of the cluster A and is related to σ 1 through 2 p − 1 = σ 1 . The negative sign of k 2 -terms implies that the diffusion always works at stabilizing the intercluster modes; the most unstable mode is the one with the smallest (but finite) wavenumber. Taking the large L limit, this smallest wavenumber is vanishingly small, so that the k 2 -terms can be dropped from the expression of the eigenvalue.
Thus the diffusion does not affect the stability, while the stability depends on the population ratio. To obtain the analytical expression of the population ratio, let us consider the interface dynamics. The two clusters A and B are treated as the fixed boundaries of the interface. Since the inter-cluster mode is stable, the boundary conditions of the interface profile are given by φ(−∞) = 0 and φ(∞) = π , and σ 2 is replaced by the steady-state value, σ 2 = 1. Then equation (4) becomes where we have defined σ ≡ σ 1 and δ ≡ β/α, and rescaled time and space as K t → t and Suppose that there is a traveling solution φ = f (x − ct), where the interface velocity c is written as c = L 2 dσ dt , owing to the fact that interface motion varies the population ratio. Multiplying (7) by ∂ x f and integrating over the entire space yields Therefore, σ has a stable solution formally written as From (8) it is seen that the interface slowly moves with the timescale of O(1/L) toward this stable state. Thus the interface dynamics, or the time evolution of σ , are decoupled from the rest. An explicit expression of the steady-state solution of σ can be obtained through perturbation expansion. First we look for a stationary solution of (7). When θ satisfies tan θ = −δ, there exists an exact solution connecting φ(−∞) = 0 and φ(∞) = π: where κ = √ 2γ cos θ . It is easily verified from (9) that this solution satisfies σ = 0. Then perturbation expansion can be performed in terms of σ up to the first order. We obtain (see appendix C) where χ(δ) = 2δ(1+δ 2 ) 1+4δ 2 coth π δ. Substituting the expression of σ into the intra-cluster eigenvalues, we obtain the stability condition of 2-cluster states. In the large L limit, the k 2 -term in the eigenvalues vanishes and the stability boundary is given by σ ± 2γ cos θ = 0.
(12) Figure 3 shows the dependence of σ on θ. We have plotted only the real part of σ , while in our simulation the imaginary part is O(K ) and is negligible. Both equations (11) and (12) fit well with the numerical data. Moreover, substituting (11) into (12) yields which gives the threshold value of θ for the stability of stationary 2-cluster states. Note that the stability condition is independent of γ . The theoretical lines given by (13) are in excellent agreement with numerical data in figure 1. Thus, within the range of parameter θ determined from (13), we can control σ as a function of θ via (11). Now the interpretation of the recurrent 2-cluster states is given as follows: given an initial condition, the system converges to a 2-cluster state with a slowly moving interface. When σ , varied through the interface motion, exceeds the threshold given by equation (12), the intracluster mode becomes unstable and one of the clusters collapses. Since the inter-cluster mode remains stable, the system returns to a 2-cluster state with reduced σ , the whole process repeated ad infinitum.

Desynchronized state
We give a theoretical analysis of equation (7) for the desynchronized state. The perfect desynchronized state is defined such that all the order parameters vanish. However, in practice, some order parameters remain finite because of the boundary effect for flux-free boundary conditions. First, as an ideal case, let us assume a periodic boundary condition. Then the perfect desynchronized state can be given by φ(x) = 2π x/L. Linear stability analysis for this profile shows that for each mode with the wavenumber k l = 2πl/L (l 1) the corresponding eigenvalue is λ (l) desync = λ (l) 0 − αk 2 l − 2βik l k 1 , where λ (1) 0 = − 1 2 , λ (2) 0 = γ e iθ and λ (l 3) 0 = 0. Hence the l = 2 mode loses stability at θ = ± π 2 in the large L limit, which fits well with the numerical data in figure 1. Note that if the diffusive coupling is absent, only l = 1 and 2 modes are stable, l 3 modes being neutral.
Since the boundary is not periodic but flux free in the present case, the profile deviates from the linear one, as shown in figure 1(c). Accordingly, the steady-state values of σ l are shifted from zero by O(1/L) (order of the width of the boundaries). Note that in the linear regime the main contribution to σ l comes from the mode k l . The modes l = 1 and 2 have the eigenvalues of order O(1) as seen above, and thus σ 1 and σ 2 remain as O(1/L). On the other hand, since l 3 modes have eigenvalues of O(1/L 2 ), nonlinear effects of order O(1/L 2 ) coming from terms such as σ 1 σ 2 make l 3 modes grow up to O(1). Therefore, in order to get a better desynchronized state, we need to add as many higher harmonics as possible, as demonstrated in section 4.

Multi-cluster states
The arguments for 2-cluster states can be extended to n-clusters in the following way. Consider the following coupling function: This coupling function, when introduced to discrete oscillators, creates stable equally populated n-clusters with the phases evenly separated (appendix B). Let us find a stationary, equally populated n-cluster solution of (4) with (14). Such a solution satisfies σ m = 0 (m < n) and σ n = 1, and hence only the nth harmonic remains in (4). Then, by choosing θ so as to satisfy we have a solution with each cluster separated by 2π n and all the n − 1 interfaces having the same interface profile given by φ(x) = 4 n arctan exp(κ n x) with κ n = √ nγ cos θ. This state is stable against interand intra-cluster fluctuations (see appendix B).
For stability of the desynchronized state, the same argument as in the above n = 2 case holds and the stability boundary is given by θ = ±π/2.

Numerical confirmation with the Brusselator model
The above analytical expressions are used for the control of oscillatory media. As a model system of oscillatory media, we adopt the Brusselator model: The As a feedback function we propose the following: where φ(u, v) is the phase of the limit-cycle oscillation 4 , and the parameters k n and ψ n are the feedback intensity and the phase shift of the nth feedback term, respectively. The coupling function is obtained from h(φ) and the phase response function Z (φ), which characterizes the sensitivity of the phase to a weak external perturbation (see appendix A). By expanding Z (φ) = l z l cos(nφ + χ l ), the coupling function is written as  (18) producing the coupling functions given by (14) with n = 2 and 5.
For l > n, k l and ψ l are equal to zero.
For n = 2 For n = 5 In principle, as long as z l is finite, we can assign any value to the lth harmonics of the coupling function by choosing appropriate values for k l and ψ l . The advantage of using (18) is that the relation between parameters in the coupling function and the feedback parameters k n , ψ n is given in a simple manner. (We could also use as the feedback h(u, v) a polynomial of u with multiple time delays [19,20], but in that case the relation is represented as a nonlinear function and the parameters need to be calculated numerically.) Hence, given a coupling function, we can calculate the corresponding feedback parameters by measuring the phase response function Z (φ). Table 1 shows Z (φ) and the feedback parameters corresponding to (14) for n = 2 and 5. In the following numerical investigation, we set γ = 0.3. First we study 2-cluster states in the one-dimensional case with the feedback parameters corresponding to n = 2. We have confirmed that, for several parameter values of θ we can observe stationary 2-clusters, recurrent 2-clusters and desynchronized states, with the order parameter values predicted by the phase model (deviation of order O(10 −2 )). As an example, in figure 3(a), numerically obtained critical values of the real part of σ 1 (denoted by '+') are superimposed on the data from the phase model, which are in good agreement with the corresponding phase model, with deviations O(10 −2 ).
Next, we use (14) for n = 5 to produce the equally populated 5-cluster state and the desynchronized state in the two-dimensional case. The 5-cluster is shown in figure 4(a), the parameter θ is given by (15) with n = 5. The order parameters are |σ 5 | = 0.739 and |σ l | ∼ O(10 −2 ) for l < 5, indicating that the five clusters are well defined and approximately equally populated. In figure 4(b) we have a desynchronized state with θ just below the threshold (θ = −1.58), where |σ l | ∼ O(10 −3 ) for l 5 and O(10 −2 ) for l > 5. Note that the degree of desynchronization becomes better than that for n = 2 shown in section 3; we can make a better desynchronized state by adding appropriate higher harmonics.

Discussion: experimental realization of the theory
To apply our method to experimental systems, we need to find the constants α and β, and the response function Z (φ). Since it is generally expected that the target pattern appears in the oscillatory media (due to inhomogeneities, or by applying a manual stimulus) [21], β can be measured by using the target pattern, assuming its phase profile to be φ(x, t) = t + k|x| (the origin is at the center of the target pattern) with the measurable quantities k and : substituting this expression into (4) with K = 0 we obtain β = ( − ω)/k 2 . The decay rate of the local perturbation from a uniform oscillation gives α. The response function Z (φ) can be measured by perturbing the system through a global parameter, to which we also apply the global feedback.
To use (18), instantaneous measurement of the phase φ(x, t) at each spatial point is needed. If at least one quantity of an oscillator is observable, this can be done by, for example, constructing a delayed coordinate.
Moreover, the following things should be taken into account for experimental realization of our theory. Firstly, the feedback must be weak for the precision of the phase description. This implies that the system size should be large enough to obtain well-defined cluster states even under weak feedback: the width of the interface is O( √ D/K ), which must be sufficiently smaller than the linear dimension L. Also, to make a coupling function containing large enough higher harmonics with weak feedback, it is better for the oscillation to be of relaxation type: then the response function has higher harmonics with large amplitudes, so that we can keep the feedback signal weak to realize a desired coupling function (see the expression for (φ) (19)). Secondly, the emergence of phase singularity leads to the breakdown of the phase description and must be avoided.
We have checked in our preliminary numerical simulations that the multi-cluster states and the desynchronized state are robust against noise. Thus we are convinced that our proposed method works in experimental systems.

Concluding remarks
We have proposed a theoretical framework for designing spatial patterns in oscillatory media. When a certain pattern is found in a phase model with a specific coupling function, the same pattern can be realized in oscillatory media by applying a properly constructed nonlinear feedback. In this paper, we found analytically tractable coupling functions that enable us to quantitatively control the spatial patterns. Using these coupling functions, we investigated the phase equation with global coupling and found the parameter regions where the following patterns stably exist: 2-cluster states with specified population ratios, equally populated multicluster states and the desynchronized state. In the case of 2-clusters, we gave an analytical expression of the population ratio of the two clusters as a function of a feedback parameter. We also proposed a simple form of the nonlinear feedback function to make the calculation of the feedback parameters easier. We exemplified all these results using the Brusselator model and succeeded in reproducing the patterns predicted by the phase model. Since our method is based on measurable quantities only, it is expected that the method will be verified in a real experiment.
The desynchronized state deserves further comment. Our results show that even in oscillatory media one can drive the system into the desynchronized state, as well as in discrete oscillators [19]. Such a control is not only of medical [19,23], but also potentially of industrial interest; for example, it would be beneficial when constant output from an oscillatory catalytic reaction is desirable.
Further investigation of the phase model with other coupling functions will be of great interest for controlling more complex patterns, although our method is limited to an oscillatory system and cannot be applied to some typical spatial patterns such as the Turing pattern. Also, investigating the control of the Benjamin-Feir unstable systems by replacing the phase diffusion equation with the Kuramoto-Sivashinsky equation will be interesting both in a theoretical sense and for applications.
In this appendix we derive the phase model (4) from a reaction-diffusion system with a global feedback. The system is assumed to undergo spatially uniform oscillation when external control is absent (namely, the system is Benjamin-Feir stable [21]). Dynamical evolution of a d-dimensional oscillatory medium is described by a reaction-diffusion equation: Note that here we consider a general situation, where the global feedback is introduced through a global parameter q. In equation (3) and in [19,20], the feedback is simply applied additively. External feedback is applied to q as where q 0 and K > 0 are constants. By assumption, ∂ t u = F(u; q 0 ) yields a limit-cycle oscillation, with its solution denoted by u = u 0 (t). The function p(t) describes a global feedback signal, given by where h(u) is some feedback function. The integration is taken over the entire space and S is the volume of the system. (Various functions can be considered for h. Our particular choice has been given in equation (18).) As we have assumed, feedback intensity K is small, so that by dropping O(K 2 ) equation (A.1) can be approximated by where f (u) ≡ (∂ F/∂q) q=q 0 . When f is independent of u, the global parameter q appears additively and the system reduces to equation (3). When a spatial pattern emerges for small K > 0, the spatial variation, and thus ∇ 2 u, is expected to be small, vanishing as K → 0. Thus, in addition to the feedback term, we may treat the diffusion term as small perturbations to the limit cycle (this is the case in our simulation, where the interface width is O( √ D/K )), and therefore the diffusion term is of the same order as the feedback). Then, following a standard method developed by Kuramoto [21], we can derive a closed description for the phase variable for our oscillatory medium. As is usually adopted, the phase φ(u) is defined so as to satisfy ∂ u φ · F(u; q 0 ) = ω. Substituting this relation into the identity ∂ t φ = ∂ u φ · ∂ t u, we obtain

B.1. Steady-state 2-cluster solution
Consider a set of N identical oscillators with the frequency ω. The dynamics are written aṡ To produce n-cluster states, it is sufficient that the coupling function contains up to the nth harmonics [22]; linear stability analysis shows that the harmonics smaller than n do not contribute to the stability of n-cluster states, and the nth harmonics works in a similar way to n : 1 periodic forcing [1]. In particular, to observe 2-cluster states, one needs to prepare the coupling function such that the first harmonics destabilizes the 1-cluster, i.e. perfect synchronization, and the second one assures the 2-cluster. Thus the coupling function for 2-cluster states can be written as The amplitude of the first harmonics can be absorbed into the coupling constant K . Also, for later convenience we choose the negative sign for the second harmonics. Assume that the oscillators form a 2-cluster state, where N A oscillators belong to the cluster A with the phase φ A and N − N A oscillators belong to the cluster B with φ B . In the phaselocking state (φ i = for all i), we get When we choose θ 1 = 0, (B.5) has a solution ψ = π for any value of p. In particular, when |γ | < 1 2 , ψ = π is the only solution except for ψ = 0, the single cluster solution. If θ 1 = 0, the phase difference is shifted from π except for p = 1 2 .

B.2. Linear stability analysis
We perform the linear stability analysis by expanding φ j as φ j = φ (0) j + ξ j , where φ (0) j = 0 for j ∈ A and φ (0) j = π for j ∈ B. First we consider the inter-cluster mode, where the fluctuation is uniform in each cluster. In this case we can write ξ j∈A = ξ A and ξ j∈B = ξ B . The mode ξ A obeyṡ Changing the variables as ξ ± = ξ A ± ξ B , and using (π) = (−π) = −1 − 2γ cos θ , we obtaiṅ The zero mode ξ + represents uniform rotation along with the limit cycle. On the other hand, ξ − corresponds to the inter-cluster mode, which is stable regardless of θ for |γ | < 1 2 . Next we consider the intra-cluster mode, where the fluctuation occurs within each cluster and the spatial average of ξ within each cluster is zero. We geṫ ξ j∈A = K (2 p − 1 − 2γ cos θ)ξ j∈A , (B.10) Thus for the intra-cluster fluctuation the eigenvalues λ A = K (2 p − 1 − 2γ cos θ) and λ B = K (1 − 2 p − 2γ cos θ ). These modes can be destabilized depending on the population ratio p.
In particular, when |θ | > π/2, either λ A or λ B is positive for any p.
Similarly, the coupling function (14) produces n-cluster solutions φ i = 2πl/n (l = 0, . . . , n − 1) for arbitrary population ratio, as can be checked by direct substitution. Stability is studied analogously to the n = 2 case above, the intraand inter-cluster eigenvalues given by where p m is the fraction of the mth cluster. When the clusters are equally populated, λ (n) intra = −K nγ cos θ and is stable for |θ | < π/2. Conversely, when |θ | > π/2 at least one of the n clusters has positive λ (n) intra and the n-cluster state is no longer stable.