Finite Size Scaling in the Kuramoto Model
We investigate the scaling properties of the order parameter and the largest nonvanishing Lyapunov exponent for the fully locked state in the Kuramoto model with a finite number of oscillators. We show that, for any finite value of , both quantities scale as with the coupling strength sufficiently close to the locking threshold . We confirm numerically these predictions for oscillator frequencies evenly spaced in the interval and additionally find that the coupling range over which this scaling is valid shrinks like with as . Away from this interval, the order parameter exhibits the infinite- behavior proposed by Pazó . We argue that the crossover between the two behaviors occurs because at the locking threshold, the upper bound of the continuous part of the spectrum of the fully locked state approaches zero as increases. Our results clarify the convergence to the limit in the Kuramoto model.
The coupled oscillator model introduced by Kuramoto in the late 70’s has established itself as a paradigmatic model for the study of synchronization phenomena, where it opened a vast area of research. The Kuramoto model allows to investigate the interplay between the tendency that individual oscillators have to run at their natural frequency and a sinusoidal all-to-all coupling which attempts to synchronize the oscillators [?]. Kuramoto elegantly solved the model in the limit of infinitely many oscillators with natural frequencies drawn from a Lorentzian distribution, for which he showed that upon increasing the coupling between oscillators, the system undergoes a transition from an incoherent disordered phase to a partially synchronized state with a finite fraction of oscillators rotating in unison [?].
Several evolutions of the original Kuramoto model have been investigated, including models with different natural frequency distributions, with couplings defined on a complex network topology, oscillators with inertia, couplings with frustration, with time delays and even negative couplings to name but a few [?]. Such extensions are motivated by the connection that the Kuramoto model has with several physical systems, ranging from synchronization phenomena in biological systems [?] to Josephson junction arrays [?], via synchronous AC electric power systems [?].
Recently, there has been a renewed interest in the finite size behavior of the Kuramoto model [?]. The problem is of interest, because all physically relevant systems and numerical simulations deal with a finite number of oscillators, which makes finding solutions to the Kuramoto model mathematically more involved. In particular, at finite , the continuum limit breaks down so that self consistent equations for physically relevant quantities such as the order parameter can no longer be written in a mathematically convenient integral form. Important steps forward in the understanding of the finite size Kuramoto model include the description of the Lyapunov spectrum for the fully locked state as well as estimates of the critical coupling necessary for synchronization to occur [?]. A complete understanding of the transition to the infinite- behavior is however still lacking.
In this work we consider the finite size Kuramoto model  on the complete graph with nodes and edges
where and are the phases and the natural frequencies of the oscillators, respectively, and is the coupling strength. For natural frequencies defined on a bounded interval, there exists a critical value of the coupling for which the system is in a fully locked state where all oscillators synchronize, with as [?]. For the particular case of uniformly distributed frequencies, the main focus of this work, it has been found that the transition from the incoherent state to full synchrony is of first order [?]. We investigate the scaling properties of the Lyapunov spectrum characterizing the linear stability of the fully locked state, and of the order parameter introduced by Kuramoto [?]. We show that above the locking threshold the largest non vanishing Lyapunov exponent scales like . Relating the expression for the order parameter [Eq. (Equation 2) below] to the Lyapunov exponents, we show that the order parameter also scales as , being the order parameter at the locking threshold. We confirm numerically these results for uniformly distributed oscillator frequencies. At first glance, our results disagree with Pazó who obtained [?]. The two results can be reconciled once one realizes that Pazó’s calculation is strictly valid for an infinite number of oscillators only, while our results are derived for finite . We find numerically that our finite result is always valid close enough to . However, its range of validity becomes narrower and narrower as increases, with numerical data consistent with . We further argue that the crossover is triggered by the dependence on of the next largest non vanishing Lyapunov exponent at , . Corrections to our results being of order , they can no longer be neglected as . A side result of our approach is that all Lyapunov exponents of the fully locked state of the Kuramoto model are monotonically decreasing functions of the coupling strength. This directly implies that the linear stability of the fully locked state improves as the oscillator coupling is increased and that if the locked state exists at , it exists at all coupling strengths . We note that this result could have been anticipated starting from the properties of the Lyapunov spectrum discussed in Ref. [?].
This paper is organized as follows. Section 2 recalls the definition of fully locked states in the Kuramoto model. Sections Section 3 and Section 4 present the calculation of the monotonicity of the Lyapunov exponents as a function of the coupling constant. Section 5 discusses the behavior of the largest nonvanishing Lyapunov exponent and of the order parameter in the immediate vicinity of the phase-locking threshold for a large but finite number of oscillators.
2The Kuramoto model
We consider the Kuramoto model defined by Eq. (Equation 1) and though our results remain valid for distributions defined on bounded intervals. Introducing the order parameter [?]
Eq. (Equation 1) can be rewritten as
Given the invariance of the Kuramoto model under a global shift of all phases, we can set . Without loss of generality we consider natural frequencies such that , which is tantamount to considering the system in a rotating frame. For , Eq. (Equation 1) admits stationary solutions given by
They are referred to as fully locked states. The linear stability of fully locked states is governed by the spectrum of the stability matrix , obtained by linearizing Eq. (Equation 1) close to , and defined as
Since the stationary solutions of Eq. (Equation 1) are invariant under a global rotation of all angles, one of the eigenvalues of is identical to zero. A stationary solution is linearly stable as long as is negative semidefinite. This condition ensures that for any small perturbation around , the system’s state, subject to the dynamics of Eq. (Equation 1), returns to exponentially fast. The eigenvalues of are referred to as the Lyapunov exponents and the linear stability condition is expressed as
In what follows , is the othonormal basis of eigenvectors of defined by . In particular is the eigenvector associated with .
According to Sylvester’s criterion, a necessary condition for to be negative semidefinite is that all its diagonal elements are negative (i.e. for all ). This implies 
The positivity of the cosine, together with Eq. (Equation 3), allows to rewrite
This choice actually corresponds to the unique stable locked state solution of the all-to-all Kuramoto model [?].
3Monotonicity of the order parameter
In this section we show that for the stable fully locked state the magnitude of the order parameter grows monotonically as the coupling constant increases. This result has already been reported in the literature , however our calculation below is based on a novel formalism which we will use later on. We therefore present it.
We start by expressing the square of the modulus of the order parameter as
where the sum runs over all pairs of oscillators. Taking the derivative of Eq. (Equation 7) with respect to gives
To obtain an expression for , we take the derivative of the stationary condition
with respect to . This gives
where and . Since is singular, we invert Eq. (Equation 9) using the Moore-Penrose pseudoinverse of defined as
where and . Multiplying Eq. (Equation 9) by yields
Finally, the difference between any two components of the expression above is given by
In order to determine the sign of the right-hand side of Eq. (Equation 13) it is useful to introduce the incidence matrix of the network. Given a graph of nodes and edges and given an arbitrary orientation of each edge, the incidence matrix is defined as follows
The product is a vector in whose entry is equal to where and are the nodes connected by edge , and where the sign of this difference depends on the arbitrary choice of orientation of the edge ( is the source and is the sink in this case). Similarly, given a vector , the product is a vector in whose entry is equal to the sum over all edges connected to node , and with the sign fixed by the nature (sink or source) of site .
We then rewrite the Kuramoto model, Eq. (Equation 1), in vector form using the incidence matrix we just introduced
where we defined for . Thus, for a stationary solution we have
This compact formulation allows to write
Injecting this last identity into Eq. (Equation 13) gives for the fully locked state
since and for in the stable fully locked state. The order parameter is therefore a monotonously increasing function of .
4Monotonicity of the Lyapunov exponents
Next, given the stable fully locked state , we compute the variation of its Lyapunov exponents as a function of the coupling strength . Because is real and symmetric, we can apply the Hellmann-Feynmann theorem [?] to calculate . We obtain
We express the derivative of the stability matrix with respect to as , which injected back into Eq. (Equation 17) yields
The matrix is defined by
Next we show that for the linearly stable fully locked state, for all values of , i.e. the Lyapunov exponents are monotonically decreasing functions of the coupling strength. If the stationary solution considered is linearly stable (i.e. ), the first term in the right-hand side of Eq. (Equation 18) is negative and only the sign of the second term needs to be determined. We note that shares the same zero row sum property as , thus and one readily obtains that as should be.
as well as
Hence, the product
is positive, since as shown in Section 3. This result implies that for the Kuramoto model on the complete graph, increasing the coupling strength systematically reduces the difference for all pairs of oscillators and .
Putting all this together, is a zero row sum matrix and Eq. (Equation 21) proves that all its off diagonal entries are positive. Thus, invoking Gershgorin’s circle theorem [?], we conclude that is negative semidefinite and thus . This concludes the proof that the Lyapunov exponents of the fully locked solution of the Kuramoto model are decreasing functions of the coupling, i.e.
This result implies that if the fully locked state is stable at , then this solution remains linearly stable and thus can be continuously followed in the interval . In other words, starting from a stable configuration, for all ensures that no instability occurs as the coupling increases (i.e. that none of the Lyapunov exponents, except , vanishes). Equation (Equation 22) not only implies that the stable fully locked state remains stable as the coupling strength is increased, but also that it becomes “more“ stable, in the sense that more negative Lyapunov exponents correspond to shorter timescales to return to equilibrium. We note that the monotonicity of the Lyapunov exponents with the coupling can also be derived starting from the properties of the spectrum of the fully locked state presented in Ref. .
5Scaling behavior of the order parameter
It is known that for uniformly distributed oscillator frequencies, the transition between the incoherent and the fully synchronized state is of first order [?]. For finite , this transition occurs as the coupling is increased above and is characterized by a discontinuous jump in the order parameter from to . The values of and depend explicitly on the number of oscillators and can be calculated for specific distributions of natural frequencies [?].
The investigation of fully locked states in the infinite- version of the Kuramoto model with frequency distributions supported on a bounded interval dates back to Ermentrout  who showed that for uniform frequency distributions the locking threshold and the order parameter at the locking transition are given by and . More recently, Pazó  showed that for the infinite- Kuramoto model and a uniform box distribution of natural frequencies, the order parameter above the locking threshold scales like
We next show analytically that for the finite size Kuramoto model and uniform frequency distribution, the scaling of the order parameter instead goes like . Furthermore, we find numerically that the range of validity of this behavior decreases with .
In Section 4, we showed that the Lyapunov exponents of the stable fully locked state decrease as the coupling is reduced. Upon reducing , the locking threshold is eventually reached at which point the locked state becomes unstable and ceases to exist. This bifurcation is accompanied with the vanishing of the Lyapunov exponent . Assuming that does not vanish at , sufficiently close to , we can approximate Eq. (Equation 23) for by
Solving Eq. (Equation 25) yields
This result indicates that the largest non zero Lyapunov exponent approaches zero with a square root behavior in the vicinity of the bifurcation. For symmetrically distributed natural frequencies , it follows from Eq. (8) in Ref. [?], together with Eq. (Equation 20) that is finite. Numerical results to be presented below for that case corroborate Eq. (Equation 26).
Eq. (Equation 23) also captures the asymptotic behavior of the Lyapunov exponents in the limit . Since the Lyapunov exponents are decreasing functions of the coupling, at large values of we have for all . Neglecting the second term in the right-hand side of Eq. (Equation 23) yields
as expected. Since decreases with for all , when the value of all cosines entering the definition of the stability matrix Eq. (Equation 4) approaches in which case its eigenvalues are with multiplicity , and with multiplicity .
We next turn our attention to the order parameter close but above locking. When the coupling approaches the locking threshold, . This justifies to truncate the sum in Eq. (Equation 16), keeping only the dominant term
Using the scaling behavior of derived above, Eq. (Equation 26), we obtain the leading expression for by replacing and respectively by and in the right-hand side of Eq. (Equation 27). Solving the resulting ordinary differential equation we obtain
To check our main results, Eqs. (Equation 26) and (Equation 28), we numerically simulate Kuramoto models with box distributed natural frequencies and various . We follow Refs. [?] and take natural frequencies evenly spaced in the interval according either to the mid-point
or the end-point rule
because they allow to obtain leading-order estimates for and , Eqs. (Equation 31) and (Equation 32) below. Few results we obtained with randomly but homogeneously distributed corroborate the results to be presented. Figs. Figure 1 and Figure 2 show numerical results for and oscillators. The data confirm the scaling predictions of Eqs. (Equation 26) and (Equation 28) sufficiently close to . Thus, we report a discrepancy between the scaling of in the thermodynamic limit and our finite size scaling which goes like . Some distance away from , one seems to recover the behavior , as is evident for .
The above reasoning predicts that the square root behavior is valid for sufficiently close to , but how close? This is investigated in Figs. Figure 3 and Figure 4 which show that the coupling range inside which the finite size scaling holds decreases with . The apparent discrepancy between Pazó’s [?] and our results is therefore the trademark of a crossover from finite to . Figure 4 gives the coupling range over which the numerical data obtained for deviates from our theoretical prediction, Eq. (Equation 28), by more than or , as a function of the inverse of the oscillator number. The observed behavior suggests that the coupling range inside which decreases with as with for both mid-point and end-point frequency distributions.
While we are not able to derive analytically the value of the exponent , we can pinpoint the origin of the crossover from to as . In our treatment above we neglected terms with in the sum over in Eq. (Equation 23). This is an increasingly bad approximation as the number of oscillators tends to infinity, because then as . To show this, we recall the finite size asymptotics recently derived in Refs.  for the mid-point and end-point frequency distributions, Eqs. (Equation 29) and (Equation 30). In Ref.  Ottino-Löffler and Strogatz obtained the finite- corrections (including numerical prefactors) of the locking thresholds
where is the Hurwitz zeta function evaluated at , and is defined by with . One obtains the leading finite size corrections to the order parameter as
Despite the different scalings for the locking threshold, the asymptotic scaling of the order parameter, Eq. (Equation 32), is the same for both frequency distributions.
Mirollo and Strogatz  further showed that the spectrum of the locked state for the finite size Kuramoto model is composed of a discrete part consisting of the eigenvalues and and of a continuous part containing the remaining eigenvalues where . An additional result of Ref.  is that for symmetric frequency distributions (as is the case for the mid-point and end-point rules) the Lyapunov exponents and can be located even more sharply as
with the second largest frequency.
At the locking threshold, and , the fully locked state is marginally stable and expanding the bounds of Eq. (Equation 33) in powers of using Eqs. (Equation 31) and (Equation 32) gives the scaling of the gap which separates the continuous part of the spectrum from zero. Using , and , respectively for mid-point and end-point rules we obtain
for both choices. Figure 5 confirms numerically the scalings of and at .
Eq. (Equation 34) then shows that at the locking threshold, . This implies that the larger the number of oscillators, the closer to zero the continuous part of the spectrum will be. Thus neglecting terms with in the sums in Eqs. (Equation 16) and (Equation 23) is an increasingly unjustified approximation as increases. The exponent in the behaviors of and relies on this truncation, which is justified only in an interval which is shrinking with . To recover the exponent obtained by Pazó in the continuous limit would require a resummation of all terms in Eqs. (Equation 16) and (Equation 23) as tends to infinity, which we have not been able to do.
We investigated the scaling properties of the Kuramoto model with uniformly distributed natural frequencies close to the synchronization threshold at finite but growing number of oscillators. We found a non trivial behavior in that both the largest non zero Lyapunov exponent , and the order parameter of the fully locked state scale like and , above the locking threshold . Our results differ from the prediction of Pazó [?] for infinitely many oscillators. We showed that this apparent disagreement is the trademark of a crossover form finite to . The range of validity of our result shrinks with . We found numerically , with . Although the numerics presented in this work are for evenly spaced frequencies, our results remain valid for other choices of ’s compatible with a uniform distribution. Our scaling predictions for and do not depend on this choice and we checked numerically on few examples that they remain valid for frequencies drawn randomly from a uniform distribution.
The fully locked states of the Kuramoto model have been thoroughly investigated in the limit of infinitely many oscillators. For the special case of uniform frequency distributions, long-established analytical results are known for: i) the value of the locking threshold, ii) the value of the order parameter at phase locking and iii) the scaling behavior of the order parameter above . For finite , however, much less is known. Finite size corrections to the locking threshold and to the order parameter for frequencies uniformly distributed over the interval have been calculated only recently [?]. The motivation behind the present work is to investigate further the finite behavior of the Kuramoto close to the locking threshold. Our manuscript complements Ref. [?] by investigating the scalings of the largest Lyapunov exponent and of the order parameter above . The observed crossover from finite to infinite and shrinking range of validity of our results significantly clarifies the mechanism behind the convergence to the limit of the Kuramoto model.
We thank S. Strogatz for useful comments. This work was supported by the Swiss National Science Foundation.