Solvable model of the collective motion of heterogeneous particles interacting on a sphere

I propose a model of mutually interacting particles on an M-dimensional unit sphere. I derive the dynamics of the particles by extending the dynamics of the Kuramoto-Sakaguchi model. The dynamics include a natural-frequency matrix, which determines the motion of a particle with no external force, and an external force vector. The position (state variable) of a particle at a given time is obtained by the projection transformation of the initial position of the particle. The same projection transformation gives the position of the particles with the same natural-frequency matrix. I show that the motion of the centre of mass of an infinite number of heterogeneous particles whose natural-frequency matrices are obtained from a class of multivariate Lorentz distribution is given by an M-dimensional ordinary differential equation in closed form. This result is an extension of the Ott-Antonsen theory.


Introduction
Systems consisting of many elements, such as a society, a school of fish, a flock of birds and the network of neurons in the brain, have attracted much research attention over years [1,2,3]. Oscillation is a ubiquitous phenomenon found in systems such as these that consist of a large number of mutually interacting elements. Other examples of these systems include reaction-diffusion systems [4], electrochemical reactions [5], electronic circuits [6,7], spiking neurons [8], the human gait [9], flashing of fireflies [10], the female menstrual cycle [11], a linear array of CO 2 waveguide lasers [12] and Josephson-junction arrays [13]. Dynamics of oscillatory element i can be described bẏ where x i is the state vector and χ i (x i ) is the function defining the behaviour of element i. The dynamics of the system consisting of these elements can be described bẏ where ψ ij (x i , x j ) is the function characterising the interaction from j to i and N is the number of oscillators in the system. In many cases, these systems are so difficult to solve that one can only perform numerical simulations and observe the results. To gain insight into the collective phenomena of these complicated systems, we need to develop methods to analyse these systems. The phase description or phase reduction of limit-cycle oscillators introduced by Kuramoto [14] is the most widely used method to analyse the synchronisation phenomena that lead to observed oscillations. The phase description describes the state of a limit-cycle oscillator by a variable φ (called the phase) whose dynamics arė where ω is the natural frequency, p(t) is an external force and Z(φ) is the phasesensitivity function. If there is no external force, this system oscillates with constant frequency ω. The external force p(t) advances or delays the phase. Z(t) determines the sensitivity of the oscillator to the external force. If two oscillators with phases φ 1 and φ 2 are coupled to each other, the force exerted by oscillator 2 on oscillator 1 is given by the phase φ 2 , i.e. p(φ 2 ). Assuming the coupling is weak, the long-time average of the mutual coupling Z(φ 1 ) · p(φ 2 ) can be regarded as a function of the phase difference φ 2 − φ 1 . Thus, the dynamics of two oscillators are described bẏ where Γ ij is the coupling function characterising the interaction from j to i.
The Kuramoto-Sakaguchi model is the simplest model of weakly coupled oscillators and is described bẏ where K is the coupling strength ( figure 1 (a)). In the Kuramoto-Sakaguchi model, the coupling function is identical among all oscillator pairs. In such a system, the natural frequency ω i differs from oscillator to oscillator; the most intensively studied is the system in which the natural frequencies follow a Lorentz distribution. The system has a critical coupling strength K c at which the synchronisation transition occurs. If K < K c , the system remains desynchronised and the order parameter r (i.e. the centre of mass) defined by is zero. For K > K c , the system is synchronised and has non-zero r. This model has been studied intensively and has contributed to the understanding of the synchronisation phenomena in general [15]. Recently, two powerful methods to analyse the Kuramoto-Sakaguchi model have attracted attention. First, Watanabe and Strogatz [16] showed that the trajectory of a system of N-coupled oscillators is on three-dimensional invariant manifold if the system is described byφ where ω is the natural frequency common to all oscillators; K, J, α, and β are time-dependent parameters. The phase of a single oscillator in the system at time t is determined by its initial phase and a nonlinear transformation involving three parameters. An important point is that the transformation with the same parameters can be used to obtain the phase of any oscillator in the system. Marvel and collaborators showed that the Möbius transformation describes the i-th phase at time t [17,18]. This transformation is a one-variable projection transformation in the complex plane where z i (t) = exp[iφ i (t)] and both of σ and a are time-dependent and oscillatorindependent parameters. This type of Möbius transformation maps the unit circle on the complex plane to itself. Because the Möbius transformation (2) has a real parameter σ and a complex parameter a, all the trajectories of the system are confined to manifolds with three dimensions. In other words, the dynamics of the system can be described by the dynamics of σ and a. This method simplifies the analysis of a system of oscillators with identical natural frequencies and enables a precise analysis of a wider range of synchronisation phenomena such as chimera states [19] and persistent fluctuations in synchronisation rates [20] than before.
Second, consider Kuramoto-Sakaguchi oscillators whose natural frequencies come from a Lorentz distribution where γ is the scale parameter. Ott and Antonsen showed that the order parameter r of an infinite number of such oscillators obeys the Stuart-Landau equatioṅ under the condition that the initial phases are on a specified manifold [21]. Ott-Antonsen theory has been used to describe systems with bimodal natural-frequency distributions [22] and systems with two subpopulations [23]. In these studies, the systems contain an infinite number of heterogeneous oscillators, but are described by differential equations with a few dynamical variables, allowing the systems to be analysed. Ott-Antonsen theory made the analysis of the synchronisation phenomena of complicated systems much simpler and much more thorough than previous methods. The Kuramoto model with α = 0 and ω = 0 can be regarded as the XY spin model if noise is injected into the oscillators. The Hamiltonian of the XY spin model where the summation runs over all connected i and j. Spins in the XY spin model are two-dimensional unit vectors. The XY spin model is generalised to the n-vector model, in which spins are unit vectors of arbitrary dimension. The XY spin model and n-vector model have been used to analyse frustrated magnets. The n-vector models are highly simplified and abstract models for describing magnetic spins; they do not necessarily correspond to physical substances. However, solving the n-vector models has advanced our understanding of phase transition. Moreover, the n-vector model with n = 0 was suggested to correspond to self-avoiding walks [24,25]. Similarly, although percolation models in high-dimensional lattices and random graphs do not necessarily correspond to physical materials, the results of infinite-dimensional percolation were useful to interpret the results of percolation on complex networks such as the attack tolerance of networks [26]. The generalisation of such models allows us to predict what happens in novel problems and to find relationships between two seemingly unrelated physical systems.
The n-vector model, which is a high-dimensional extension of the XY spin model, has contributed to statistical physics by improving our understanding of phase transitions. A high-dimensional extension of the Kuramoto model may similarly advance the understanding of the collective phenomena of coupled systems and can be applied to other fields. If the extended model is solvable, it would have many possible applications. In particular, if the methods similar to those of Watanabe and Strogatz and Ott and Antonsen are applicable to the model, it would be quite useful for understanding collective phenomena.
The study of collective phenomena goes back a long time. In the 17th century, Huygens observed the antiphase locking of two pendulum clocks [27]; in the 18th century, Kaempfer reported the synchronised flashing of fireflies at the banks of the Chao Phraya [10]. Although these descriptions date back to the early modern period, the study of synchronisation phenomena advanced only in the latter half of the 20th century with the development of the phase-description method of limit-cycle oscillators [14,28]. Similarly, a solvable generalisation of the Kuramoto-Sakaguchi model would be useful for understanding the collective phenomena in a broader setting.
Motivated by this historical perspective, I propose an extension of the Kuramoto-Sakaguchi model to describe the collective motion of particles interacting on a unit sphere ( figure 1 (b)). First, by extending the framework of the Watanabe-Strogatz transformation, we derive the dynamics of the individual elements, which I call 'particles'. The present states of these particles, which I sometimes call 'positions', are provided by the projection transformation (a high-dimensional extension of Möbius transformation) of their initial values. This dynamics of a particle are given bẏ where † denotes the Hermitian conjugate, x is an M-dimensional real or complex vector representing the state of the particle, Ω is an M × M matrix corresponding to the natural frequency of phase oscillators and g determines the force exerted on the particle. Assuming that x is a two-dimensional real vector and assuming that x is a one-dimensional complex vector lead to the same dynamics as the Kuramoto-Sakaguchi model. The position of the particle at time t is given by the same projection transformation, irrespective of the initial position x 0 . Second, we derive the centre of mass r of the particles that initially are uniformly distributed on the sphere. The centre of mass corresponds to the order parameter of the Kuramoto-Sakaguchi model. In particular, I show that the centre of mass is given by one of the parameter vectors of the projection transformation if the state vector x is a complex vector. Third, we extend the Ott-Antonsen theory to high-dimensional systems. By assuming that the naturalfrequency matrices of particles are obtained from a multivariate Lorentz distribution, I show that the centre of mass of particles with complex variables is described by lowdimensional ordinary differential equations. In particular, I show that, if the particles in a system are attracted to Kr, where K is an M × M matrix, the dynamics of r are described byṙ whereΩ is a matrix determined by the probability distribution of the natural frequencies of the particles. We derive the limit cycle from these dynamics and show that a transition emerges from a desynchronised state to a synchronised state that is similar to that found in the Kuramoto model. I also show that this theoretically derived limit cycle agrees well with the results of numerical simulation. Finally, I discuss the possible applications of the present model and problems to be addressed in the future.

Dynamics of the particles on the unit sphere induced by the projection transformation
The Möbius transformation (2) has been shown to underlie the dynamics of (1) [18]. Marvel and collaborators showed that the phase of oscillator i at time t is given by the Möbius transformation (2) of its phase at time 0. Notably, given the initial conditions of oscillators, the same parameters σ and a (i.e. the same Möbius transformation) can be used to calculate the phases at t of all oscillators in the system. To extend the Kuramoto-Sakaguchi model to an M-dimensional system, we use the Möbius transformation. Because this transformation is a one-dimensional projection transformation, we derive the dynamics of variables whose present values are given by an M-dimensional projection transformation of the initial values. We consider a system in which the M-dimensional state vector x of a particle is given by the projection transformation of the initial state x 0 , where the M × M matrix A, M-dimensional vectors b and c and scalar value d are time dependent. The system can be real valued or complex valued (i.e. x can be a real or a complex vector). To derive the dynamics of particles whose time evolution is described by (3), we differentiate x with respect to time to obtaiṅ Substituting into Qx 0 + q, we have Replacing Q and q with c † and d, respectively, yields The first term in the numerator of (4) is and the second term in the numerator of (4) is Rearranging terms giveṡ The dynamicsẋ can be obtained from (5) if To realise (6), the dynamics of the parameters of the projection transformation must bė Thus, there are two ways to obtain the state x at time t. First, to directly calculate (6) from the initial condition x 0 . Second, to calculate the parameters of the transformation, A, b, c and d, using (7a)-(7d) and the initial condition A = I, b = 0, c = 0 and d = 1 and then map x 0 to x using (3). Note that, with these initial conditions for the parameters, A, b, c and d, the relationship x = x 0 is satisfied at t = 0. The first and the second method calculate the time evolution of the state vector and the transformation, respectively. Thus, one can regard these methods as corresponding to the Schödinger and Heisenberg representation of quantum mechanics. The dynamics of phase oscillators can be regarded as those restricted to a unit circle in a plane. Similarly, we restrict the dynamics of a particle to the sphere |x| = 1. The condition is satisfied by the constraints Thus, the dynamics are described bẏ where g is an arbitrary vector and Ω is an antisymmetric matrix in the real-valued system and an anti-Hermitian matrix in the complex-valued system. This is a subclass of the differential equations called the Riccati matrix differential equation [29,17]. For a real two-dimensional vector reduces toẋ where we set and Replacing x 1 with cos φ and x 2 with sin φ gives where we assume g 1 = K cos ψ and g 2 = K sin φ, which is the dynamics of a Kuramoto oscillator under an external force. Similarly, for a complex variable x = [exp(iφ)], (9) reduces toφ where Im denotes the imaginary part and we set Ω = (iω) and g = [g] = [K exp(iψ)/2], which is the dynamics of a Kuramoto oscillator as well. Thus, these two systems have the same dynamics. However, in general, a system of 2M-dimensional real state variables and a system of M-dimensional complex state variables do not have the same dynamics because the antisymmetric matrix Ω has 2M 2 − M real degrees of freedom in the former case, whereas the anti-Hermitian matrix Ω has M 2 real degrees of freedom in the latter case. These two examples suggest that g and Ω correspond to the external force and the natural frequency, respectively, in the Kuramoto-Sakaguchi model. The naturalfrequency matrix Ω determines the direction and the speed of the rotation of the particle on the unit sphere. Next, we consider whether the interaction terms −xg † x + g of (9) can be regarded as the tangential component of a central force on the surface of the unit sphere. Let us assume that the central force is described by whose origin is g. Assuming that x † x = 1, the tangential component is calculated by using the projection operator I − xx † as which gives −xg † x + g if x and g are real vectors and f (|g − x|) = |g − x|. Thus, the interaction of particles in this system is the same as that for the system of objects that are attracted to g by ideal springs and are restricted to lie on the M-dimensional unit sphere. Potential energy is obtained by integrating the force −xg † x + g. Assuming that particle i is attracted to particle j (i.e. g = x j ) the force exerted on particle i by particle j is given by −x i x † j x i + x j . The potential energy is given by where we notice that the position vector x on the sphere is orthogonal to the tangent vector dx. E is the potential energy of the n-vector model.

Dynamics of the parameters of the projection transformation
Constrained by (8a) and (8b), the dynamics of the parameters arė where we set d = 1. In this system, the parameters A, b and c of the projection transformation are dependent variables. The projection transformation to and from the unit sphere satisfies From (11a), c is determined by A and b. Substituting (11a) into (11b), we obtain Rearranging the terms and multiplying by A † −1 from the left and by A −1 from the right, we have

The matrix A satisfying this equation is defined by
where U is an arbitrary orthogonal matrix if x is a real vector and an arbitrary unitary matrix if x is a complex vector and with Σ 1/2 = ( Σ ij ) and V ΣV † being the singular value decomposition of H = (1 − |b| 2 )I + bb † .
V appears two times in the singular value decomposition of H because H is a Hermitian matrix. Taken together, these relationships imply that the projection transformation to and from the unit sphere is determined by a vector b and an orthogonal or unitary matrix U . Because the degrees of freedom of an M-dimensional orthogonal matrix and an M-dimensional unitary matrix are M(M − 1)/2 and M 2 , respectively, the degree of freedom of the transformation is M(M + 1)/2 in the real-valued system and M 2 + 2M in the complex-valued system. Thus, the real-valued systems with M = 2 have three degrees of freedom, and the complex-valued systems with M = 1 also have three degrees of freedom. It is identical to the dimension of the invariant manifolds of Watanabe and Strogatz [16]. The vector b is an eigenvector of H with eigenvalue 1 because Vectors orthogonal to b (i.e. vectors n satisfying b † n = 0) are eigenvectors of H with eigenvalues 1 − |b| 2 because Hn = (1 − |b| 2 )n + bb † n = (1 − |b| 2 )n.
H 1/2 approaches I in the limit |b| → 0, which indicates that the initial conditions A = I and b = 0 are consistent with each other if we set U = I.

Centre of mass of particles
We set the initial conditions of the parameters of the transformation to A = I and b = 0 to satisfy x = x 0 . The projection transformation of x 0 to x is given as We now derive the centre of mass m = 1 N 1≤i≤N x i of the particle. We assume that, initially, the points are uniformly distributed on the unit sphere. This assumption simplifies the calculation of the centre of mass and enables the low-dimensional description of the entire system. Assuming that the points are uniformly distributed on the unit sphere in the initial conditions and that there are an infinite number of particles, we have for real-valued systems and x dx 0 for complex-valued systems, where S M = 2π M/2 Γ(M/2) is the surface area of the M-dimensional unit sphere. Equation (14) suggests that the density of points and the centre of mass are completely determined by b because U does not affect the distribution after transformation if the particles are uniformly distributed on the unit sphere. Thus, without loss of generality, we can replace U of (14) with I for calculating the centre of mass m. Decomposing x 0 into x 0 , and n is a vector orthogonal to b (i.e. n = x 0 − ηb 1 ) we have where we use (12) and (13) and set U = I. The centre of mass of the particles in real space is given by where 2 F 1 (a, b; c; x) is the ordinary hypergeometric function. The centre of mass of the particles in complex space for M = 1 is given by because n = 0 in this case. For M ≥ 2, the centre of mass is given by Thus, b is the centre of mass of the system with complex variables. Figure 2 shows the results of simulation for particles described by a real threedimensional variable. In this system, all particles have the natural-frequency matrix

Results of simulation for the system with identical particles
and experience the external force No mutual coupling among particles is introduced in figure 2(a). The solid line in figure 2(a) is obtained by the direct simulation of (9) and shows the trajectory of one of the particles. The dashed line shows the trajectory obtained by simulating the dynamics of parameters A and b by (10a) and (10b) and transforming the initial conditions of the particle by the projection transformation (3) with these parameters. These two simulation methods yield the same result. Adding a mutual interaction and setting the external force to where m is the centre of mass of N = 1000 particles and I perform the direct simulation with N = 1000 particles, which initially are randomly distributed on the unit sphere. The same natural-frequency matrix as (16) is used for all particles. The solid line in figure 2(b) shows the trajectory of the centre of mass of the particles in the system. Approximating m of (17) by (15) (Figure 2(b), dashed line), I simulate the dynamics of parameters A and b by (10a) and (10b) and obtain the trajectory with the projection transformation (3). Figure 2(b) shows that the direct simulation of the original system agrees well with this approximation. This agreement occurs because the initial distribution of particles may be regarded as uniform when the number of particles is sufficiently large.

Extension of the results of Ott and Antonsen
Next, we consider a heterogeneous system in which the natural-frequency matrix Ω varies from particle to particle. Ott and Antonsen showed that the order parameter (i.e. the centre of mass) of an infinite number of Kuramoto-Sakaguchi oscillators obeys the Stuart-Landau equation if their natural frequencies are obtained from a Lorentz distribution [21]. To extend their results, we introduce two additional assumptions. First, we assume a complex-valued system and that the state x of the particle is a complex vector. Because the centre of mass m is approximated by b in a complexvalued system, this assumption facilitates the derivation of the dynamics of the centre of mass. Second, we assume that the natural-frequency matrix Ω is drawn from a multivariate Lorentz distribution. We assume that the entries of the natural-frequency matrix Ω are given by the linear superposition of random variables drawn from Lorentz distributions: where ζ z is obtained from the Lorentz distribution with location parameter µ z and scale parameter γ z > 0, whose probability density is given by Because Ω must be an anti-Hermitian matrix, the eigenvalues of Ω z (0 ≤ z ≤ Z) are purely imaginary. We assume that where K(t) is the coupling strength, f (t) is the external force, and p(ζ) = 1≤z≤Z p z (ζ z ). Here, we show that, if Ω 0 is an anti-Hermitian matrix and the eigenvalues of Ω z are positive and imaginary for z > 0, the dynamics of this heterogeneous system are described bẏ where First, we derive the analytic continuation of m(Ω, t), assuming that ζ z -s are complex numbers. Because m can be replaced by b in complex-valued systems, we obtain from (10b) the dynamics of m(Ω, t), which iṡ m(Ω, t) = −m(Ω, t)g(t) † m(Ω, t) + Ωm(Ω, t) + g(t).
Assuming that m(Ω, t) satisfies the Cauchy-Riemann equations where ζ z = u z + iv z and u z and v z are real numbers, thenṁ(Ω, t) satisfies the Cauchy- where we used These relations indicate that m(Ω, t) is an analytic function of ζ z if m(Ω, t ′ ) is an analytic function of ζ z , where t ′ < t [30]. Here, we have Because m(Ω, t) is an analytic function, we have where C is a semicircle in the upper half of the complex plane with radius S and centred at the origin, and if the second integral on the right-hand side converges to zero. Because Ω z is an anti-Hermitian matrix, the time evolution of |m(Ω, t)| 2 on |m(Ω, t)| = 1 satisfies where Re denotes the real part. If all eigenvalues of the Hermitian matrix iΩ z are negative on the sphere, d dt |m(Ω, t)| 2 ≤ 0, then |m(Ω, t)| = 1 when ζ k is in the upper-half complex plane (i.e. v z ≥ 0). Thus, m(Ω, t) remains finite if Ω z are anti-Hermitian matrices with positive imaginary eigenvalues and Im ζ k ≥ 0. In the limit of large S, the dynamics of m(Ω, t), where |ζ k | = S, can be approximated bẏ where we approximate Ω by ζ k Ω k and ignore the terms without Ω. If we use the same assumption that all eigenvalues of Ω k are positive and imaginary, then the real part of all eigenvalues of ζ k Ω k are negative under the condition Im ζ k > 0. In this case, m(Ω, t) approaches to zero as S increases, and so the second integral of (20) converges to zero in the limit of large S. Thus, we have Because m(Ω k , t) is an analytic function of ζ z (z = k), the integration can be done recursively to obtain (19), so the dynamics of the order parameter r(t) are given by (18). Note that |m(Ω , t)| remains finite and converges to zero in the limit of large S = |ζ z | (z = k) even if Ω is replaced byΩ k .

Limit-cycle oscillation of the centre of mass
Assuming that there is no external force (i.e. f (t) = 0) and that the mutual interactions among oscillators are constant (i.e. K(t) = K), the dynamics are described bẏ which is a high-dimensional extension of the Stuart-Landau equation. The stability of the desynchronised state, r(t) = 0, is determined by the eigenvalues ofΩ + K.
Because Ω 0 is an anti-Hermitian matrix, Ω z (z > 0) are anti-Hermitian matrices with positive imaginary eigenvalues and γ z > 0, so the real part of is non-positive. Therefore, all real parts of eigenvectors ofΩ are non-positive, so the system without any mutual interaction (K = 0) remains desynchronised. IfΩ + K has an eigenvector e i with eigenvalue λ i , where Re λ i > 0, the limit cycle is given by and |e i | = 1. Because we get The real part of the right-hand side is positive because Re λ i > 0. Re(e † iΩ † e i ) ≤ 0 because the real part of x †Ω x is non-positive for any x. Thus, we have Re(e † i K † e i ) > 0, and so (22) is satisfied with a positive R 2 i and a λ i with a positive real component. Assuming that λ 1 is the eigenvalue with the largest real component, we prove that only the limit cycle r 1 is stable. Let us examine the time evolution of the perturbed solution where d = 1≤j≤M, j =i d j e j , |d i | ≪ 1, |ρ| ≪ 1 and |θ| ≪ 1. Because we have replacing d with 1≤j≤M, j =i d j e j , we obtain where we used R 2 i e † i K † e i = λ i − iξ i and (Ω + K)e i = λe i . Thus, the dynamics of d j is given byḋ Re λ 1 > Re λ i . Hence, the limit cycle r i is unstable if i > 1. Conversely, d j for j > 1 is stable if i = 1. This fact allows us to set d = 0 in examining the stability of the limit cycle r 1 . The resulting dynamics of ρ and θρ Hence, ρ has a stable fixed-point at ρ = 0, and θ is neutrally stable.
The solid lines are the trajectories of the order parameter r of N = 10 000 particles obtained from the direct simulation oḟ where the order parameter is calculated by The dashed lines show the trajectories of r obtained from the reduced dynamics (21). The results from the reduced dynamics agree quite well with those of the direct simulation. By varying the value of k, I obtain the fixed point and the limit-cycle oscillation of r ( figure 3(a) and (b)). This simple, low-dimensional behaviour of the limit cycle is in sharp contrast with the complicated motions of individual particles. Figure 3(c) shows the trajectories of two particles in the system of figure 3(b). One of the particles is entrained into the collective synchronisation whereas the other particle has a complicated trajectory. Figure 3(d) shows how the radius of the limit cycle depends on k. The radius obtained by solving (22) (dashed line) agrees well with that obtained by the direct simulation (crosses) and by the reduced dynamics (circles).

Conclusions
In this study, we have extended the Kuramoto-Sakaguchi model to model particles interacting on a high-dimensional unit sphere. The dynamics are described by a type of matrix Riccati differential equation, which is characterised by an external force vector and a natural-frequency matrix. The position of a particle at a given time is obtained by the projection transformation of the initial position. This result is an extension of the Watanabe-Strogatz theory. By assuming that the particles are uniformly distributed on a unit sphere in the initial conditions, the centre of mass of the particles with the same natural-frequency matrix is determined by the vector of the parameters of the projection transformation. In a system of particles with complex variables and naturalfrequency matrices with a multivariate Lorentz distribution, the motion of the centre of mass of all the particles can be described by a high-dimensional extension of the Stuart-Landau equation. This result is an extension of the Ott-Antonsen theory. A periodic solution of the extended Stuart-Landau equation agrees with the motion of the centre of mass of the system. Thus, we have shown that the collective motion of a system composed of elements with large degrees of freedom can be reduced to the dynamics of a low-dimensional system. Although the present model is an extension to high degrees of freedom of the Kuramoto-Sakaguchi model, it differs from the previously proposed extensions [31,32]. For example, Ritort introduced a variable into the model to keep the oscillators on the unit sphere [31], whereas the particles of the present model (described by (9)) remain on the unit sphere without additional terms. In addition, unlike the dynamic variables of the model by Gu and coworkers, which are represented by matrices [32], the dynamic variables in the present model are represented by vectors. Finally, unlike the present model with heterogeneous particles, neither of these two models has been reported to be reducible to low-dimensional systems.
The present model provides a method to study new types of collective phenomena by reducing the system behaviour of elements with large degrees of freedom to that of a low-dimensional system. The present model and its low-dimensional description can be applied to problems already studied from the viewpoint of Kuramoto oscillators [15], such as time-delay systems [33,34], systems of multiple-peak natural frequencies [22], non-local coupling [35], dynamics on complex networks [36], associative memory [37] and common-input synchronisation [38]. In particular, because the system of oscillators with heterogeneous interaction delays [33] and the system of oscillators whose natural frequencies obey a mixture of Lorentz distributions [22] have been studied by using the Ott-Antonsen theory, I expect that the high-dimensional extension of these systems can be solved by the present method.
I have not found a physical system whose dynamics are described by the present model. However, the phenomena observed in the present model might be useful in interpreting experimental observations of the systems that can be regarded as a population of particles on a sphere. In particular, the normalised velocity of birds in a flock approximated by the Heisenberg model [3] could be analysed by using the real-valued systems with M = 3. The results of the model suggest that the state transition from the desynchronised state to synchronised state can occur in systems with heterogeneous particles ( figure 3(d)). The present results also suggest that the trajectory of the centre of mass can be a very simple limit cycle (figure 3(b)) even if individual particles in the system exhibit complicated trajectories ( figure 3(c)). The present model could be used to approximate and to analyse systems exhibiting these properties.