Center manifold reduction for large populations of globally coupled phase oscillators

Center manifold reduction for large populations of globally coupled phase oscillators

Abstract

A bifurcation theory for a system of globally coupled phase oscillators is developed based on the theory of rigged Hilbert spaces. It is shown that there exists a finite-dimensional center manifold on a space of generalized functions. The dynamics on the manifold is derived for any coupling functions. When the coupling function is , a bifurcation diagram conjectured by Kuramoto is rigorously obtained. When it is not , a new type of bifurcation phenomenon is found due to the discontinuity of the projection operator to the center subspace.

pacs:
05.45.Xt

Introduction. Collective synchronization phenomena are observed in a variety of areas, such as chemical reactions, engineering circuits, and biological populations Pikovsky et al. (2001). In order to investigate such phenomena, a system of globally coupled phase oscillators of the following form is often used:

 dθidt=ωi+KNN∑j=1f(θj−θi),i=1,⋯,N, (1)

where denotes the phase of an -th oscillator, denotes its natural frequency drawn from some distribution function , is the coupling strength, and is a -periodic function. When , it is referred to as the Kuramoto model Kuramoto (1984). In this case, it is numerically observed that if is sufficiently large, some or all of the oscillators tend to rotate at the same velocity on average, which is referred to as synchronization Pikovsky et al. (2001); Strogatz (2000). In order to evaluate whether synchronization occurs, Kuramoto introduced the order parameter , which is given by

 r(t)e√−1ψ(t):=1NN∑j=1e√−1θj(t). (2)

When a synchronous state is formed, takes a positive value. Indeed, based on some formal calculations, Kuramoto assumed a bifurcation diagram of : Suppose . If is an even and unimodal function such that , then the bifurcation diagram of is as in Fig.1(a). In other words, if the coupling strength is smaller than , then is asymptotically stable. If exceeds , then a stable synchronous state emerges. Near the transition point , is of order . See Strogatz (2000) for Kuramoto’s discussion.

In the last two decades, several studies have been performed in an attempt to confirm Kuramoto’s assumption. Daido Daido (1996) calculated steady states of Eq. (1) for any using an argument similar to Kuramoto’s. Although he obtained various bifurcation diagrams, the stability of solutions was not demonstrated. In order to investigate the stability of steady states, Strogatz and Mirollo and coworker Mirollo and Strogatz (1990, 2007); Strogatz and Mirollo (1991); Strogatz et al. (1992) performed a linearized analysis. The linear operator , which is obtained by linearizing the Kuramoto model, has a continuous spectrum on the imaginary axis. Nevertheless, they found that the steady states can be asymptotically stable because of the existence of resonance poles on the left half plane Strogatz et al. (1992). Since the results of Strogatz and Mirollo and coworker are based on the linearized analysis, the effects of nonlinear terms are neglected. However, investigating nonlinear bifurcations is more difficult because has a continuous spectrum on the imaginary axis, that is, a center manifold in the usual sense is of infinite dimension. In order to avoid this difficulty, Crawford and Davies Crawford and Davies (1999) added a noise of strength to the Kuramoto model. The continuous spectrum then moves to the left side by , and thus the usual center manifold reduction is applicable. However, their method is not valid when . An eigenfunction of associated with a center subspace diverges as because an eigenvalue on the imaginary axis is embedded in the continuous spectrum as . Recently, Ott and Antonsen Ott and Antonsen (2008) found a special solution of the Kuramoto model, which allows the dimension of the system to be reduced. Their method is applicable only for because their method relies on a certain symmetry of the system Marvel et al. (2009). Furthermore, the reduced system is still of infinite dimension, except for the case in which is a rational function. Thus, a unified bifurcation theory for globally coupled phase oscillators is required.

In the present letter, a correct center manifold reduction is proposed by means of the theory of generalized functions, which is applicable for any coupling function . It is shown that there exists a finite-dimensional center manifold on a space of generalized functions, despite the fact that the continuous spectrum lies on the imaginary axis. This will be demonstrated for two cases, (I) and (II) , and two distribution functions , (a) a Gaussian distribution and (b) a rational function (e.g., Lorentzian distribution ). For (I), we rigorously prove Kuramoto’s conjecture, and for (II), a different bifurcation diagram will be obtained, as was predicted by Daido Daido (1996). The different bifurcation structure is shown to be caused by the discontinuity of the projection to the generalized center subspace. All omitted proofs are given in Chiba ().

Settings. The continuous model of Eq. (1) is given by with , where is defined to be

 ηl(t)=∫R∫2π0e√−1lθρt(θ,ω)g(ω)dθdω,

and is a probability measure on parameterized by . In particular, is a continuous version of Kuramoto’s order parameter. Setting yields

 dZjdt = √−1jωZj+√−1jKfjηj+√−1jK∑l≠jflηlZj−l (3) := TjZj+√−1jK∑l≠jflηlZj−l, ηj(t) = ∫RZj(t,ω)g(ω)dω=(¯¯¯¯¯¯Zj,P0)=(P0,Zj), (4)

where is the inner product on the weighted Lebesgue space , and . The linear operator is known to have a continuous spectrum on the imaginary axis. Furthermore, there exists a positive constant , such that if , has eigenvalues on the right half plane (such that the de-synchronous state is unstable), while if , has no eigenvalues. For example, if is an odd function and if is even and unimodal, then . In the present letter, for simplicity, we assume that (for , which is true if and only if ). When , has no eigenvalues, and thus the dynamics of the linearized system is quite nontrivial. In Chiba (), the spectral theory on rigged Hilbert spaces is developed to reveal the dynamics of the linearized system.

A rigged Hilbert space consists of three spaces : a space of test functions, a Hilbert space , and the dual space of (a space of continuous anti-linear functionals on , each element of which is referred to as a generalized function). We use Dirac’s notation, where for and , is denoted by . For , we have . The canonical inclusion is defined as follows. For , we denote by , which is defined to be

 i(ψ)(ϕ)=⟨ϕ|ψ⟩:=(ϕ,ψ)=∫R¯¯¯¯¯¯¯¯¯¯¯ϕ(ω)ψ(ω)g(ω)dω. (5)

By the canonical inclusion, Eq. (3) is rewritten as an evolution equation on the dual space as

 ddt|Zj⟩=T×j|Zj⟩+√−1jK∑l≠jfl⟨P0|Zl⟩⋅|Zj−l⟩, (6)

where is a dual operator of defined through .

Here, the strategy for the bifurcation theory of globally coupled phase oscillators is to use the space of generalized functions rather than a space of usual functions. The reason for this is explained intuitively as follows. If we use the space to investigate the dynamics, then the behavior of itself will be obtained. However, it is neutrally stable because of the conservation law . What we would like to know is the dynamics of the moments of , in particular, the order parameter. This suggests that we should use a different topology for the stability of . (Note that the definition of stability depends on definition of the topology.) For the purpose of the present study, is said to be convergent to as if and only if

 ∫R∫2π0ϕ(ω)e√−1jθdρt(θ,ω)→∫R∫2π0ϕ(ω)e√−1jθd^ρ(θ,ω)

for any and . The topology induced by this convergence is referred to as the weak topology. By the completion of with respect to the weak topology, we obtain a space of generalized functions . On the space , a function converges as if and only if converges for any . Since the order parameter is written as , this topology is suitable for the purpose of the present study. For this topology, it turns out that is not neutrally stable.

A suitable choice of depends on . When is a Gaussian distribution, let be the set of holomorphic functions defined near the upper half plane such that is finite. Set . We can introduce a suitable topology on so that the dual space becomes a complete metric space, which allows the existence of a center manifold on to be proven. When is a rational function, is a space of bounded holomorphic functions on the real axis and the upper half plane. In this case, we can show that is a finite-dimensional vector space, which implies that Eq. (6) is essentially a finite-dimensional system. This is why in Ott and Antonsen (2008); Martens et al. (2009), the system is reduced to a finite-dimension system when is the (sum of the) Lorentzian distribution.

Let be the semigroup generated by . In Chiba (), the spectral decomposition of is obtained by means of the rigged Hilbert space. Define resonance poles of to be roots of the equation

 ∫R1λ−√−1ωg(ω)dω+2πg(−√−1λ)=2K, (7)

on the imaginary axis and the left half plane. Roughly speaking, a resonance pole is a continuation of an eigenvalue to the second Riemann sheet of the resolvent  Strogatz et al. (1992); Chiba (). In the following, for simplicity, we assume that all ’s are single roots of Eq. (7). Define a functional to be

 ⟨ϕ|μj⟩=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∫R¯¯¯¯¯¯¯¯¯¯¯ϕ(ω)g(ω)λj−√−1ωdω+2π¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕ(−√−1λj)g(−√−1λj),(Re(λj)<0),limx→+0∫R¯¯¯¯¯¯¯¯¯¯¯ϕ(ω)g(ω)(x+√−1yj)−√−1ωdω,(λj=√−1yj).

The ’s are called the generalized eigenfunctions associated with the resonance poles due to the equality . Then, we can prove that the semigroup is expressed as

 (eT1t)×|ψ⟩=∞∑n=0K2Dneλnt⟨ψ|μn⟩⋅|μn⟩, (8)

for any , which gives the spectral decomposition of on the dual space , where are constants defined by

 Dn=limλ→λn1λ−λn(1−K2∫Rg(ω)λ−√−1ωdω−πKg(−√−1λ)).

(In the previous paragraph, was chosen so that the right-hand side of Eq.(8) converges.) Equation (8) completely determines the dynamics of the order parameter for the linearized system. If , then all resonance poles lie on the left half plane. As a result, decays to zero exponentially as , which proves the asymptotic stability of the de-synchronous state.

Center manifold reduction. When , there exist resonance poles on the imaginary axis. The generalized center subspace is defined as a space spanned by generalized eigenfunctions associated with resonance poles on the imaginary axis, say . Equation (8) suggests that the projection to is given by

 Πc|ψ⟩=M∑n=0K2Dn⟨ψ|μn⟩⋅|μn⟩. (9)

In general, is continuous only on a subspace of because the topology on is too weak. When , it is proven in Chiba () that is continuous only on . For a solution of Eq. (3), are included in , although are not. Because of the discontinuity of , an interesting bifurcation occurs when . In what follows, assume that is an even and unimodal function. In this case, on the imaginary axis, there is only one resonance pole at . Hence, is of one dimension, where is the generalized eigenfunction associated with . Next, let us derive the dynamics on a center manifold. The derivation is performed in the same way for both (a) a Gaussian distribution and (b) rational functions.

(I) First, we consider . In this case, equations for do not depend on . Thus, acts continuously on solutions of Eq. (6). Set . Then, Eq. (6) for is rewritten as

 ddt|Z1⟩=T×10|Z1⟩+ε2⟨P0|Z1⟩|P0⟩−K2⟨P0|Z1⟩|Z2⟩, (10)

where is defined by replacing with in the definition of . In order to obtain the dynamics on the center manifold, using the spectral decomposition, we set

 |Z1⟩=Kc2α|μ0⟩+|Y1⟩,α(t)=1D0⟨Z1|μ0⟩, (11)

along the direct sum . The purpose here is to derive the dynamics of . Since and are included in the stable subspace, the center manifold theorem Chiba () reveals that on the center manifold, . In particular, the last term of Eq. (10) is of order . Apply the projection to both sides of Eq. (10). Noting that is also of order because the projection is continuous on , we obtain the dynamics on the center manifold as

 ddtα=αD0Kc(ε+πg′′(0)K4c16|α|2)+O(εα2,ε2α,ε3,α4). (12)

If , this equation has a fixed point of order when . It is easy to verify that the order parameter is given by Thus, the dynamics of the order parameter is also given by Eq. (12). Since , when is even and unimodal, (de-synchronous state) is unstable, and the nontrivial fixed point (synchronous state) is asymptotically stable when , which confirms Kuramoto’s diagram.

(II) Assume that with . Then, Eq. (6) for is given by

 ddt|Z1⟩=T×10|Z1⟩+ε2⟨P0|Z1⟩|P0⟩ (13) − K2(⟨P0|Z1⟩|Z2⟩+h⟨P0|Z2⟩|Z3⟩−h⟨P0|Z2⟩|Z−1⟩).

In this case, , on which is discontinuous, appears. As before, satisfies and since acts continuously on . This implies that the term in Eq. (13) yields a cubic nonlinearity as case (I). On the other hand, since is discontinuous on , we can show that even if . As a result, the last term in Eq. (13) yields a quadratic nonlinearity. Indeed, we can verify that

 Π0|Z−1⟩=πKcg(0)D0e−√−1arg(α)+O(α). (14)

Applying the projection to both sides of Eq. (13), we obtain the dynamics on the center manifold as

 dαdt=αD0Kc(ε−K3cCh2(1−h)αe−√−1arg(α))+h.o.t., (15)

where is a negative constant, which proves that there exists a fixed point that is expressed as

 |η1|∼|α|=2(1−h)K3cCh(K−Kc)+O((K−Kc)2). (16)

Therefore, for , a stable branch emerges when , and for , an unstable branch emerges when (Fig.1(b)).

Discussion. Equation (6) shows that the dynamics of is independent of if and only if . In other words, Eq. (6) splits into two systems: a system of and a system of . Since the projection is continuous on a solution of the former system, we can show the existence of a smooth center manifold. Note that Eq. (1) is invariant under the rotation on a circle. As a result, the dynamics on the center manifold is also invariant under the rotation . If a center manifold is smooth, then the dynamics on this manifold with the rotation symmetry must be of the form . Thus, a cubic nonlinearity is dominant, and a pitchfork bifurcation generally occurs, as shown in Eq. (12). On the other hand, if , then the equations of depend on , on which is not continuous. In such a case, the center manifold is not smooth, and quadratic nonlinearity may appear, as described above. In this manner, different bifurcations occur when . Although the diagram shown in Fig.1(b) looks like a transcritical bifurcation, Eq. (15) is different from the normal form of a transcritical bifurcation. Because of the factor caused by the discontinuity of , Eq. (15) remains invariant under the rotation despite the existence of a quadratic nonlinearity. The discontinuity induces a new type of bifurcation including .

A center manifold reduction for globally coupled phase oscillators was also developed by Crawford and Davies Crawford and Davies (1999) with a noise of strength . Although they also expected a diagram such as shown as Fig.1(b) when , the factor was not obtained. Since the eigenfunction diverges as , expressions of the dynamics on the center manifold were not shown explicitly. In the present letter, we have shown that the eigenfunction exists on a space of generalized functions, which provides a correct center manifold reduction. The diagram shown in Fig.1(b) was also obtained by Daido Daido (1996) by means of a self-consistent analysis. Unfortunately, his results were not correct because he performed inappropriate termwise integrations of certain infinite series. According to his results, the order parameter is given as , which suggests that some degeneracy occurs when . However, the numerical results given in Fig.2 show that the critical exponent of the order parameter changes only when , which agrees with the results of the present study (16). Ott and Antonsen Ott and Antonsen (2008) found an inertia manifold given by when . The center manifold of the present study is a finite-dimensional submanifold of the inertia manifold, which provides a further reduction of the results of Ott and Antonsen. The key strategy of the present theory is to use spaces of generalized functions and the weak topology. The weak topology is suitable for investigating the dynamics of moments of probability density functions. Since the strategy is independent of the details of the models, this strategy will be extended to various types of large populations of coupled systems and evolution equations of density functions, such as the Vlasov equation.

Acknowledgements.
The present study was supported by Grant-in-Aid for Young Scientists (B), No.22740069 from MEXT Japan.

References

1. A. Pikovsky, M. Rosenblum,  and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, Cambridge, 2001).
2. Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984).
3. S. H. Strogatz, Physica D, 143, 1 (2000).
4. H. Daido, Physica D, 91, 24 (1996).
5. R. E. Mirollo and S. H. Strogatz, J. Stat. Phys., 60, 245 (1990).
6. R. E. Mirollo and S. H. Strogatz, J. Nonlinear Sci., 17, 309 (2007).
7. S. H. Strogatz and R. E. Mirollo, J. Stat. Phys., 63, 613 (1991).
8. S. H. Strogatz, R. E. Mirollo,  and P. C. Matthews, Phys. Rev. Lett., 68, 2730 (1992).
9. J. D. Crawford and K. T. R. Davies, Physica D, 125, 1 (1999).
10. E. Ott and T. M. Antonsen, Chaos, 18, 037113 (2008).
11. S. A. Marvel, R. E. Mirollo,  and S. H. Strogatz, Chaos, 19, 043104 (2009).
12. H. Chiba, “A proof of the Kuramoto’s conjecture for a bifurcation structure of the infinite dimensional Kuramoto model,” (submitted), arXiv:1008.0249.
13. E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So,  and T. M. Antonsen, Phys. Rev. E, 79, 026204 (2009).
14. R. Tönjes, N. Masuda,  and H. Kori, Chaos, 20, 033108 (2010).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters