Order parameter analysis for low-dimensional behaviors of coupled phase-oscillators

# Order parameter analysis for low-dimensional behaviors of coupled phase-oscillators

Jian Gao Department of Physics and the Beijing-Hong Kong-Singapore Joint Centre for Nonlinear and Complex Systems (Beijing), Beijing Normal University, Beijing 100875, China    Can Xu Department of Physics and the Beijing-Hong Kong-Singapore Joint Centre for Nonlinear and Complex Systems (Beijing), Beijing Normal University, Beijing 100875, China    Yuting Sun Department of Physics and the Beijing-Hong Kong-Singapore Joint Centre for Nonlinear and Complex Systems (Beijing), Beijing Normal University, Beijing 100875, China    Zhigang Zheng Department of Physics and the Beijing-Hong Kong-Singapore Joint Centre for Nonlinear and Complex Systems (Beijing), Beijing Normal University, Beijing 100875, China College of Information Science and Engineering, Huaqiao University, Xiamen 361021, China
July 15, 2019
###### Abstract

Coupled phase-oscillators are important models related to synchronization. Recently, Ott-Antonsen(OA) ansatz is developed and used to get low-dimensional collective behaviors in coupled oscillator systems. In this paper, we develop a simple and concise approach based on the equations of order parameters, namely, order parameter analysis, with which we point out that the OA ansatz is rooted in the dynamical symmetry of the order parameters. With our approach the scope of the OA ansatz is identified as two conditions, i.e., infinite size of the system and only three nonzero Fourier coefficients of the coupling function. Coinciding with each of the conditions, a distinctive system out of the scope is taken into account and discussed with the order parameter analysis. Two approximation methods are introduced respectively, namely the ensemble approach and the dominating-term assumption.

## I Introduction

Understanding the intrinsic mechanism of collective behaviors of coupled units has become a focus for a variety of fields, such as biological neurons, circadian rhythm, chemically reacting cells, and even social systems kuramoto1984chemical (); acebron2005kuramoto (); strogatz2000kuramoto (); pikovsky2002synchronization (); dorogovtsev2008critical (); arenas2008synchronization (); zheng1998phase (). Some properties of collective behaviors depend on the complexity of the system, while the other properties, such as phase transition, may be described by low-dimensional dynamics with macroscopic variables watanabe1993integrability (); watanabe1994constants (); marvel2009identical (); marvel2009invariant (); ott2008low (). Discovering the method to simplify the system is just as important and as fascinating as the discovery of the complexity of it.

Like most cases in physics, simplification and low-dimensional reduction are associated with some symmetry of the system. As the identity of gas particles is the foundation of statistical mechanics and collective variables as temperature and pressure, the identity of the coupled units in a complex system is also related to some order parameters. In previous works, in the limit of large number of oscillators with special coupling function this work has been done by Ott-Antonsen(OA) ansatz ott2008low () for oscillators with nonidentical parameters. As for the oscillators with identical parameters the same result was also got from group theory analysis in marvel2009identical () called Watanabe-Strogatz’s approach.

In this paper, we will show a different way in which the low-dimensional reduction is a natural consequence of the symmetry of the system by taking the order parameters as collective variables, namely the order parameter analysis. Our approach is simple and concise for oscillators with both identical and nonidentical parameters. We can also get the scope of the OA ansatz. Two more cases beyond the scope are further discussed with appropriate approximations. With our approach we show that the OA ansatz can be used beyond its scope with the approximations works.

The model discussed in this paper is the all-connected coupled phase-oscillators. In the first section, the dynamical equations for order parameters are derived. The OA ansatz is got naturally from the symmetry of the order parameter equations, with its scope as the limit of infinitely many oscillators and the condition that only three Fourier coefficients of the coupling function are nonzero. In the second and third sections, we will consider two approximate use of this ansatz beyond its scope, the case of a finite-size system and the case of coupled oscillator systems with more complicated coupling functions, with approximations respectively, i.e., the ensemble approach and the dominating-term assumption.

## Ii Equations for order parameters

The famous Kuramoto model for the process of synchronization attracts much attentions upon it is proposed and has been developed for decades. This model consists of a population of N coupled phase oscillators with natural frequencies , and the dynamics are governed by

 ˙φj(t)=ωj−N∑k=1Ajksin(φk−φj),   j=1,…,N. (1)

With the mean-field coupling and the definition of order parameter

 α=1NN∑j=1eiφj, (2)

the equations Eq. (1) can be rewritten as

 ˙φj(t)=ωj−K2i(αe−iφj−¯αeiφj),   j=1,…,N, (3)

where is the imaginary unit, is the complex conjugate of . Apart form parameters and , the dynamics of each phase variable depends only on itself and . It is the important character of this mean-field model, and the order parameter is always used to describe the state of the system, as the system is in synchronous states if and only if .

A more general form of the mean-field model with phase oscillators can be written as

 ˙φj(t)=F(α,φj,β,γj),  j=1,…,N, (4)

where is any smooth, real, -periodic function for . is the order parameter with the th order parameter defined as

 αn=1NN∑j=1einφj. (5)

is the identical parameter such as the coupling strength which is identical for all the oscillators. is the nonidentical parameter such as the natural-frequency which is nonidentical for the oscillators and usually has a distribution among the system.

Almost all the mean-field models based on the Kuramoto model belong to this general category. In the following, we will build our approach for this general model to explore the conditions in which we can get the low-dimensional description of the system Eq. (4).

### ii.1 Identical oscillators

To begin with, let us consider the simple case of identical oscillators, e.g., for in the Kuramoto model, which reads

 ˙φj(t)=F(α,φj,β),j=1,…,N. (6)

In the limit of infinitely many oscillators , let denote the fraction of oscillators that lie between and . Because each oscillator in Eq.(6) moves with the angular velocity , the single oscillator density obeys the continuity equation as

 ∂ρ∂t+∂(ρv)∂φ=0, (7) v=F(α,φ,β).

If the phase density is known, all of macroscopic properties of the system can be got through some statistical average, such as order parameters as

 αn=∫2π0ρ(t,φ)einφdφ,n=1,2,…. (8)

If we are only concerned with the collective or macroscopic state of the system, as in this paper, the macroscopic description of the phase density is equivalent to the microscopic description of phases of all the oscillators.

Moreover, for the coupling function is always -periodic for , we have the Fourier expansion of ,

 F(α,φ,β)=∞∑j=−∞fj(α,β)eijφ. (9)

With the the definition Eq. (8), the dynamics of order parameters can be got as

 ˙αn =∫2π0ineinφ˙φdφ (10) =∫2π0ineinφF(α,φ,β)dφ,

where . Substituting the expansion of into these equations, we have the closed equations for the dynamics of order parameters as

 ˙αn=in∞∑j=−∞fj(α,β)αj+n,n=1,2,…. (11)

On the other hand, from Eq. (8) the order parameters are exactly Fourier components of the phase density. If all the order parameters are known, we have

 ρ(t,φ)=12π∞∑j=−∞αje−ijφ. (12)

Hence if we know all the initial values of order parameters, with the dynamical equations Eq. (11) and the Fourier transformation Eq. (12), the system is identified, which is the same as the dynamical equations Eq. (6) for phase variables and the continuity equation Eq. (7) for the phase density . Therefore, we can choose either the order parameters, the phase variables or the phase density to perform the analysis of the collective behaviors of the coupled phase oscillators.

These three descriptions of the system, i.e. phase variables, density of phase and order parameters, correspond to the dynamical, statistical and macroscopic descriptions of the system, respectively, and they are equivalent to each other in the limit of infinitely many oscillators . As a matter of fact, the Ott-Antonsen ansatz ott2008low () is based on the representation of the density of phase, and the Watanabe-Strogatz’s approach marvel2009identical () is based on the scenario of the phase variables. In the following, we will take our approach on the base of order parameters.

The complexity of dynamical equations for order parameters Eq. (11) depends on the coupling function, or explicitly, on the Fourier expansion of . First, let us consider the simplest case, with only the first three nonzero terms of the Fourier expansion,

 F(α,φ)=f1(α)eiφ+f−1(α)e−iφ+f0(α). (13)

Because is a real function, we have where means the complex conjugate of . With Eq. (13), the dynamics of order parameters Eq. (11) becomes

 ˙αn=in(f1αn+1+¯f1αn−1+f0αn), (14)

with and . The recursion form of these equations shows that there is some structure for order parameters with which we can simplify the system and get some low-dimensional equations for the high-dimensional coupled oscillator dynamics of the system. In fact, by choosing an invariant manifold as , all the equations for is reduced to

 ˙α1=i(f1(α1)α21+¯f1(α1)+f0(α1)α1) (15)

for . The infinitely many dynamical equations of order parameters is thus reduced to a single equation on this manifold and the corresponding phase density defined by Eq. (12) is the so-called Poisson kernel distribution as

 ρ(t,φ)=12π1−r21−2rcos(φ−ϕ)+r2, (16)

where . If we choose the initial state on this manifold, the state will never evolve out of it, which is called the invariant manifold of the dynamical system Eq. (14). This is exactly the low-dimensional behavior of the system we are looking for. With our approach, the manifold is derived naturally and concisely. We will call this manifold the Poisson manifold in this paper as in marvel2009identical () where the dynamics of , as Eq. (15), is got from the group theory analysis for Josephson junction arrays.

### ii.2 Nonidentical parameters

Furthermore, the low-dimensional behavior is not confined to the system of identical oscillators. The order parameter analysis we used above can also be used in a more general case of nonidentical oscillators with nonidentical parameters.

Firstly, let us consider a system of nonidentical oscillators with a discrete distribution of the parameter . This distribution naturally separates the oscillators into groups, and each group has the same parameter . In this case, the phase variable for oscillators is denoted by which means the th oscillator with the same parameter .

For each group we can define a local order parameter as

 α(r)n=1nrnr∑l=1einφ(r)l,n=1,2,…, (17)

where is the number of oscillators in the group with . The order parameter of the system denoted by is defined by

 αn=p∑r=1α(r)nnrn,n=1,2,…. (18)

Assume that the coupling function for each oscillator depends only on the order parameters whether for the whole system or the groups, then equations for phase oscillators are

 ˙φ(r)j(t)=F(α,α(r),φ(r)j,ω(r),β),j=1,…,nr, (19)

where . and are the local and global order parameters, respectively.

As the coupling function is -periodic for , by performing Fourier expansion of , the equations for order parameters of each group read

 ˙α(r)n=in∞∑j=−∞fj(α,β,α(r),ω(r))α(r)j+n, (20)

where , and .

The equations Eq. (20) are closed for local order parameters, and the local order parameters here can be regarded as the coordinates of the system which are equivalent to phase variables. Each group is almost separated from others and the only relation for the groups is the dependence of the coupling function on order parameters of the system.

In the simplest case when only the first three terms of Fourier expansion are nonzero, we have

 ˙α(r)n=in(f1α(r)n+1+¯f1α(r)n−1+f0α(r)n), (21)

with and . For each group, is obviously a solution for the equations which reduce all the equations for to the same one as

 ˙α(r)1=i(f1(α(r)1)2+¯f1+f0α(r)1), (22)

where . The Poisson manifold is an invariant manifold for each group, and the order parameters of the system read

 αn=p∑r=1(α(r)1)nnrn,n=1,2,…, (23)

which is determined by . Therefore, the Fourier coefficients of the coupling function depends on only. The system described by the local order parameters is governed by a group of Poisson manifolds with the low-dimensional dynamics Eq. (22).

In the limit with and , the distribution of parameter as for has the continuous form denoted by the function and the local order parameters become the function of , i.e., . Replacing the summation by the integral over , we get the continuous form of the dynamical equations as

 ˙α1(ω)= i(f1(ω,α1(ω),α∗)(α1(ω))2+¯f1(ω,α1(ω),α∗) (24) +f0(ω,α1(ω),α∗)α1(ω)),

where are functionals of as

 α∗n=∫(α1(ω))ng(ω)dω,n=1,2,…. (25)

Even though for each specific the approach has already reduced the dynamics of the oscillators to the dynamic of as Eq. (24), it is still hard to get any analytical results as Eq. (24) depends on the integral Eq. (25) which cannot be expressed by functions of . On the other hand, for some specific choice of , the integral Eq. (25) can be obtained analytically, in which case we can get the simpler form of Eq. (24).

What we are looking for is a two-step reduction. The first step is finished by introducing the Poisson manifold with which we have reduced the dynamic of phase oscillators to the dynamic of a group of Poisson manifolds, i.e., Eq. (24). The second step is rooted in the relation between the Poisson manifolds in the group which will reduce the dynamics of the group further to the dynamic of a specific Poisson manifold in this group.

For instance, in the case that the function is analytical which can be extended to the complex plane, and the distribution of is the Lorentzain distribution as

 g(ω)=Δπ((ω−ω0)2+Δ2), (26)

the integral Eq. (25) can be obtained by residue theorem as . Setting , Eq. (24) becomes

 ˙α1(−i)= (27) +f0(−i,α1(−i))α1(−i)),

which is closed for order parameter . Then, the behavior of , together with the collective behaviors of the system, is determined by Eq. (27).

Another solvable example is the case when the distribution of natural frequencies of oscillators is the Dirac’s delta distribution function . The integral Eq. (25) can be worked out as . Setting in Eq. (24), the model is reduced to the network of identical oscillators which we discussed above.

The approach shown above indicates there are two steps of the reduction scheme. The first step is the reduction from the equation Eq. (20) to Eq. (22) or the continuous form Eq. (24), which means the choice of the Poisson manifold of the system. This is related to the the symmetry of the system or the recursion form of order parameters equations. The second step is the further operation from Eq. (24) to Eq. (27), which depends on the special choice of distribution of nonidentical parameters. This reduction with the Lorentzain distribution was first found in ott2008low () with an ansatz for low-dimensional manifold, namely the OA ansatz. For the case of the Dirac’s delta distribution function, this reduction was discussed comprehensively in marvel2009identical () with group theory analysis of the system. In our approach, the reduction comes from the same basis, i.e., the symmetry of order parameter equations.

Along the approach, we can also consider more nonidentical parameters, such as the location of oscillators for example in Laing2009chimera (), where the coupling function depends on the locations as with the distributions and . The order parameter of the system is defined as

 αn=∫∫αn(ω,x)g(ω)q(x)dωdx,n=1,2,…. (28)

When the Fourier expansion of the coupling function for contains only the first three terms, for each specific and , the relation is exactly the solution for dynamical equations, which reads

 ˙αn(ω,x)=in(f1αn+1(ω,x)+¯f1αn−1(ω,x)+f0α∗n(ω,x)). (29)

If is analytic for and is the Lorentzain distribution, the integral Eq. (28) can be obtained for as

 αn =∫∫αn(ω,x)g(ω)q(x)dωdx (30) =∫αn(−i,x)q(x)dx,

where . Setting and in the dynamical equation Eq. (29), we have

 ˙γ1(x)= i(f1(ω,γ1(x),α)γ1(x)2 (31) +¯f1(ω,γ1(x),α) +f0(ω,γ1(x),α)γ1(x)),

where is the order parameter of the system and the functional of as

 αn=∫(γ1(x))nq(x)dx,n=1,2,…. (32)

The solution of this integral differential equation Eq. (31) describes the structure of the local order parameter along .

Up to now, we have discussed the system of all-connected phase oscillators in the limit of infinitely many oscillators and the case that only three Fourier coefficients of the coupling function are nonzero. The low-dimensional invariant manifold, namely the Poisson manifold, is got for both cases. We will see in the next two sections that these two conditions are exactly the scope of the OA ansatz. Two cases beyond this cope will be discussed respectively.

## Iii Network with finite size

We have considered the general model Eq. (6) in the limit of infinitely many oscillators, , in which we could get the phase density of oscillators and the corresponding continuity equation Eq. (7). Order parameters could be considered as Fourier coefficients of , from which we get the equivalent expression of the dynamics of the system as the order parameter equations Eq. (11) and build the approach above.

In the limit , the macroscopic variable, namely the order parameter , has the limit for steady states as the synchronous state and incoherent states, which gives us the basis to discuss the system analytically. In the case of a finite but large number of oscillators, i.e., , the order parameter defined as will fluctuate around the value . This fluctuation depends on the number of oscillators as Daido1986 (); Daido1989 (); Daido1990 (). When the fluctuation is small enough, as for , the collective behaviors of the system with finite number of oscillators could be described by the approximation of infinitely many oscillators, for which analytical methods can be applied. Some analytical methods, e.g., the self-consistent method and the OA ansatz, are applicable for the case of infinitely many oscillators or equivalently the phase density description.

On the other hand, in the case of only a few number of oscillators, e.g., , the difference would be so large, as , whose magnitude comparable with the value of . Obviously, it is not reliable to treat the system with the approximation in this case. It is necessary to propose new approaches in analytically dealing with the collective behaviors of finite-oscillator systems.

### iii.1 Ensemble approach

In our approach, order parameters could be considered as not only the Fourier coefficients of but also the collective variables as

 αn=1NN∑j=1einφj, (33)

which does not depend on the approximation of infinitely many oscillators. For the system as

 ˙φj(t)=F(α,φj,β),j=1,…,N, (34)

with the definition of order parameters Eq. (33), the dynamical equations for order parameters are

 ˙αn=inKK∑j=1einφj˙φj,n=1,2,…. (35)

For the coupling function is -periodic for , we have the Fourier expansion as

 F(α,φ,β)=∞∑j=−∞fj(α,β)eijφ. (36)

Substituting the expansion into Eq. (35), together with Eq. (33), we have the closed form of equations for order parameters as

 ˙αn=in∞∑j=−∞fj(α,β)αj+n,n=1,2,…, (37)

which is the same as we get in the limit of infinitely many oscillators.

In the case of a few number of oscillators, the definition of order parameters can also be considered as the coordinate transformation that transforms the microscopic variables to macroscopic variables . Therefore it is inspiring that the dynamics of phase oscillators with corresponding initial conditions is equivalent to the dynamics of order parameters with corresponding initial conditions . The number of oscillators, no matter finite or infinite, has no influence on the dynamical equations of order parameters. This forms the important basis of our approach.

Following this approach, when only the first three Fourier coefficients of the coupling function are nonzero, by setting and substituting the relation to the dynamical equations, the equations Eq. (37) for all the will be reduced to a single one as

 ˙α1=i(f1(α1)α21+¯f1(α1)+f0(α1)α1). (38)

It appears that we could get the Poisson manifold again even the size of the system is finite. However, the finite-size effect should be taken into account with care.

As a matter of fact, in the case of a finite number of oscillators, the Poisson manifold with is not attainable for the system. To see this, take the system of as an example. For the first two order parameters and , by definition

 α1=122∑j=1eiφj,α2=122∑j=1e2iφj, (39)

if the system described by is in the Poisson manifold, the relation could be rewritten in terms of the definition Eq. (33) as

 122∑j=1e2iφj=(122∑j=1eiφj)2. (40)

However, the equality Eq. (40) is not naturally valid. By using a simple calculation of the difference between the left-hand and the right-hand terms in Eq. (40) one obtains

 122∑j=1e2iφj−(122∑j=1eiφj)2=14(eiφ1−eiφ2)2. (41)

This indicates that the system can evolve on the Poisson manifold only when . In general, the relation gives exactly infinitely many independent constraints like Eq. (40) with , and any solutions for finite oscillators will be determined as the trivial one as for all the index , which means that except the synchronous state, all the states of the system of finite oscillators lie out of this Poisson manifold. The system can only evolve on the Poisson manifold in the limit of infinitely many oscillators.

Whereas, as a matter of fact, the equations Eqs. (37) hold for all , whether the size of system is finite or infinite, and the synchronous state always satisfies the relation of the manifold. The Poisson manifold can be used at least approximately in the vicinity of synchronous state. However, we need some new methods and approximations.

Let us consider an ensemble with identical systems of oscillators which have the same dynamical equations and same parameters. The initial phases of oscillators are chosen from the same distribution, which makes the mean values of order parameters over the ensemble have the limit when . This leads to the definition of the ensemble order parameter as

 ⟨˙αn⟩=limM→∞1MM∑j=1α(j)n,n=1,2,…, (42)

where is the number of sampling systems in the ensemble and is the th order parameter for the th system in the ensemble. Taking the ensemble average for both sides of the dynamical equations Eq. (37), we get

 ⟨˙αn⟩=in(⟨f1αn+1⟩+⟨¯f1αn−1⟩+⟨f0αn⟩), (43)

where , and means ensemble average.

In general Eq. (43) is not solvable because the terms in the right side as cannot be simply described by in general. However, if is independent of all the systems in the ensemble as , then we can get exactly the dynamical equations of the ensemble order parameters, as

 d⟨αn⟩dt=in(f1⟨αn+1⟩+¯f1⟨αn−1⟩+f0⟨αn⟩), (44)

where and the terms may depend on . Similar to Eq. (37), there is a solution for the ensemble order parameter equations Eqs. (44) as , which is exactly the Poisson manifold.

Moreover, for more general cases , but when the terms can be approximated by , namely the statistical independence, we can also get the dynamical equations of the ensemble order parameters with similar approach. The validity of this approach can be measured by the error terms

 ejn=⟨fjαn⟩−⟨fj⟩⟨αn⟩, (45)

with and .

### iii.2 The star Sakaguchi-Kuramoto model

In the following, let us take the star Sakaguchi-Kuramoto model as an example, which is a typical topology and model in grasping the essential properties of heterogeneous networks and synchronization process, as

 ˙θh =ωh+λN∑j=1sin(θj−θh−β), (46) ˙θj =ω+λsin(θh−θj−β),

where , and are the natural frequency and the phase of hub and leaf nodes respectively, is the coupling strength and is the phase shift. By introducing the phase difference , the dynamical equation can be transformed into

 ˙φj=Δω−λN∑k=1sin(φk+β)−λsin(φj−β), (47)

where is the natural-frequency difference between hub and leaf nodes.

By introducing the order parameter , we could rewrite the dynamics as

 ˙φj=feiφj+g+¯fe−iφj,j=1,⋯,N, (48)

where , . Hence the system is in the framework which we discussed above, with the first three nonzero Fourier coefficients of the coupling function. With our approach, the dynamical equations for order parameters can be obtained as

 ˙αn=in(fαn+1+¯fαn−1+gαn),n=1,2,…. (49)

In this specific model, we cannot simply make the approximation as because diverges with . There are two ways in further simplifying the system. First, we could set a hypothetical system as

 ˙ϑj=Δω−NλM∑j=11Msin(ϑj+β)−λsin(ϑj−β), (50)

where , is the hypothetical phase oscillator and is considered as a parameter in this system. Defining the order parameter for this system as , Eq. (50) becomes

 ˙φj=feiφj+g+¯fe−iφj,j=1,⋯,M, (51)

where , . Following our approach, the corresponding order parameter equations read

 ˙αn=in(fαn+1+¯fαn−1+gαn),n=1,2,…, (52)

which are exactly the same as Eq. (49). But for this hypothetical system, we could take the limit of infinitely many oscillators as , which could be discussed analytically with the traditional OA anstaz.

Obviously, these two systems are quite different for phase variables, one with finite oscillators the other with infinitely many oscillators, but they have the same order parameter equations. In the following we will see that the hypothetical system Eq. (50) is exactly a representation of the ensemble of the original model Eq. (47).

On the other hand, following our ensemble approach, for the system Eq. (47), let us choose an ensemble consisting of systems with the same parameters and , and different initial conditions chosen from the same distribution, as the Poisson kernel. We have the ensemble average of the dynamical equation as

 ⟨˙αn⟩ =in(iλ2e−iβ⟨αn+1⟩−iλ2eiβ⟨αn−1⟩+Δω⟨αn⟩ (53) +i2λN(⟨α1αn⟩eiβ−⟨¯α1αn⟩e−iβ)),

where . Setting

 ⟨α1αn⟩=⟨α1⟩⟨αn⟩+en1, (54) ⟨¯α1αn⟩=⟨¯α1⟩⟨αn⟩+en2,

for . If and , or typically and , they can be seen as the perturbation terms for the dynamics of ensemble order parameters in Eq. (53). Ignoring the perturbations, we could get

 ⟨˙αn⟩ =in(iλ2e−iβ⟨αn+1⟩−iλ2eiβ⟨αn−1⟩ (55) Δω⟨αn⟩+i2λN(⟨α1⟩⟨αn⟩eiβ −⟨¯α1⟩⟨αn⟩e−iβ)),

and this is the same as Eq. (52) for hypothetical system. In this case, the ensemble order parameters are the same as the order parameters for the hypothetical system.

For the case of a finite but large number of oscillators, the hypothetical system is exactly the continuous model which shares the same order parameter equations but with infinitely many oscillators. As in the case of a few oscillators, the hypothetical system is a presentation of the ensemble for the systems, where the difference terms and by this method will not only introduce some fluctuations but also some systematic errors, which can be described by

 En =|⟨α1αn⟩−⟨α1⟩⟨αn⟩|, (56) E′n =|⟨¯α1αn⟩−⟨¯α1⟩⟨αn⟩|,

with for this specific model.

In terms of numerical computation, we can check the above approximation by examining and for the ensemble of systems. The approximation and can be divided into two parts as and where and are the fluctuation parts proportional to which depend on the size of the ensemble and and are the systematic parts which are introduced by the statistical independence assumption.

Take as an example. When the fluctuation parts are ignorable. We find that except for is around , all of and for are much smaller than in this case as plotted in Fig. 1(a) and (b). Comparing the typical magnitude of values of order parameters as , the approximation of the ensemble approach is obviously reasonable.

For the initial conditions chosen from the Poisson kernel, with the dynamical equations Eq. (55), the ensemble order parameters will evolve on an invariant manifold, i.e., the Poisson manifold with . Denote by , we have

 ˙z=−λ2z2e−iβ+i(Δω+i2λN(zeiβ−¯ze−iβ))z+λ2eiβ. (57)

The difference between and the ensemble average of numerical simulations is checked in Fig. 1(d), which shows that it is reasonable to describe the behavior of the ensemble by . Moreover, for every single system in the ensemble, with some fluctuation, it can also be described by the ensemble order parameter and hence by , as shown in Fig. 1(c).

With this ensemble assumption works, the two-dimensional equation Eq. (57) in the bounded space as describes the dynamics of the coupled oscillator system. Every stationary state for phase oscillators system has its counterpart in this space for . For example, the in-phase state(IPS) defined as corresponds to a limit cycle as . The synchronous state (SS) defined as corresponds to a fixed point with and . The splay state(SPS) defined as with the period of corresponds to a fixed point with and . In the two-dimensional order-parameter space, it is easy to analyze the existence and stability conditions of these states, and the basion of attraction for each state in the coexistence region can also be conveniently determined xu (). This is the reason why we introduce the invariant manifold governed by low-dimensional dynamics.

As we see, in the case with only finite oscillators, the ensemble order parameter plays the dominant role in revealing the key collective structure of the system. The price of the ensemble order parameter description is the error induced by the statistical independence assumption, whose validity can be checked by numerical results. In the next section we will see that when more than three Fourier coefficients of the coupling function are nonzero, where the traditional order parameter scheme fails, one can still discuss the collective dynamics in terms of the ensemble order parameter approach.

## Iv Complex coupling function

In the above sections, for the coupling functions with only first three terms as

 F(α,φ)=f1(α)eiφ+f−1(α)e−iφ+f0(α), (58)

we have the dynamical equations for order parameters as

 ˙αn=in(f1αn+1+¯f1αn−1+f0αn),n=1,2,…. (59)

On the invariant manifold , all these equations are reduced to a single one

 ˙α1=i(f1(α1)α21+¯f1(α1)+f0(α1)α1). (60)

For general situations, the coupling function can not be simply truncated to only the first terms. It is thus important to consider the order parameter approach for more complicated cases.

### iv.1 Higher order coupling

Firstly, let us consider a coupling function with higher order Fourier coefficients like

 G(α′,ϕ)=gm(α′)eimϕ+g−m(α′)e−imϕ+g0(α′), (61)

where is an integer and are the order parameters for the system. Note that Eq. (61) can be regarded as the transformation of Eq. (58) with and . The order parameters for is related to the order parameters for as

 αn=1KK∑j=1einφj=1KK∑j=1ein⋅mϕj=α′n⋅m. (62)

And the phase density of with the transformation reads

 ρ′(ϕ) =ρ(φ) (63) =12π∞∑j=−∞αjeijφ=12π∞∑j=−∞α′j⋅meij⋅mϕ.

This indicates that all the terms with are zero. Moreover, we can transform the Poisson manifold for to the invariant manifold for , which is governed by

 ˙α′m=im(f′m(α′m)2+¯f′m+f′0α′m), (64)

corresponding to Eq. (60) for . This is exactly the low-dimensional manifold for the coupling function Eq. (61) with higher order Fourier coefficients. This manifold is based on the th order parameter and the statement that all the terms with are zero. This is a bit queer in this respect.

Let us take the case as an example Skardal2011higer (). In this case only and among the Fourier components of the coupling function are nonzero, i.e.,

 F(α,φ)=f2(α)e2iφ+f−2(α)e−2iφ+f0(α). (65)

According to the above analysis, the system has the invariant manifold as , or equivalently , with dynamics

 ˙α2=2i(f2α22+¯f2+f0α2). (66)

The phase density of oscillators in this case is

 ρ(φ)=12π∞∑j=−∞α2⋅je2ijϕ. (67)

By using our order-parameter-analysis approach, for the case of , it is not hard to get the the equations of order parameters as

 ˙α2n+1 =i(2n+1)(f2α2n+3+¯f2α2n−1+f0α2n+1), (68) ˙α2n =i(2n)(f2α2n+2+¯f2α2n−2+f0α2n),

with . The manifold with is obviously the solution of these equations. This invariant manifold separates the relation into two parts, where the part conserves the relation, and the odd part can be regarded as a special choice of the relation. Even if in the state defined as , we still have , which represents the state of cluster synchrony as discussed in Skardal2011higer (). And the manifold with