Introduction

Globally coupled phase oscillators offer a benchmark description in a large number of physical and biological systems, such as synchronizing fireflies, neurons in human brain, cardiac pacemaker cells, power grids and Josephson junction arrays, just to name a few examples. In his seminal work, Kuramoto introduced a solvable model for synchronization of globally coupled phase oscillators1, which inspired and stimulated a wealth of extensive studies along the last three decades2, starting from the early investigations on the stability of its coherent and incoherent states3,4.

Almost all works on the Kuramoto model are focused on describing and characterizing the continuous phase transition (PT) which marks the passage from its incoherent state to its (partially) coherent state, i.e., a second-order PT. However, Yeung and Strogatz reported that a first-order PT with an associated hysteresis loop occurs in a Kuramoto model with time delay5. Furthermore, abrupt transitions to synchronized states of networked Kuramoto-like oscillators were also reported for a particular realization of a uniform frequency distribution (evenly spaced frequencies) and an all-to-all network topology6 and later on the same finding was described for both periodic and chaotic phase oscillators in a heterogeneous degree-distribution with positive correlations between the node degree and the corresponding oscillator's natural frequency7,8. Such a latter phenomenon was usually termed as explosive synchronization (ES) in literatures and the investigation on ES in generalized Kuramoto models has recently attracted a great attention9,10,11,12,13,14,15,16,17,18. In particular, Ref. 19 introduced a specific generalization of the Kuramoto model, in which the coupling strength is frequency-weighted.

In this paper, we provide the full analytical treatment of the model of Ref. 19 and we present a series of exact results, which include: (i) the demonstration of the first-order character of the PT for different frequency distributions (FDs), such as triangle, Lorentzian, Gaussian, bimodal Lorentzian, etc., (ii) the analytical solutions for the critical coupling strengths for both the forward and backward PT under typical FDs, (iii) the closed form of the forward transition point and the linear stability analysis for the incoherent state for Lorentzian FD and (iv) the closed forms for both the stable and unstable coherent states (and their stabilities) for the backward transition, which are rigorously demonstrated not to depend on the specific form of the FDs.

We start by introducing the generalized Kuramoto model that was the object of the study of Ref. 19:

Here θi and ωi are the instantaneous phase and the natural frequency of the ith oscillator, dot denotes a temporal derivative and κ is the coupling strength. The set of N natural frequencies {ωi} is drawn from as given FD with density g(ω). In this study, without loss of generality, g(ω) is assumed to be an even function with g(ω) = g(−ω), i.e., symmetric and centered at 0. The most important characteristic in the above model is that a weight factor ωi is assigned to the coupling. This can be understood as a weight to the uniform coupling strength κ. The introducing of this weight is significant in certain networks. For example, in social networks, different individuals may have different strength of response to the same influence/stimuli, e.g., rumors from friends. If we roughly describe the human dynamics by Kuramoto phase oscillators, a weighted factor in terms of the natural frequency of oscillators can reasonably characterize the heterogeneity of couplings in such network20. In Ref. 7, ES was first found in a special scale-free network, i.e., Barabási-Albert (BA) network, when a strictly positive correlation is required between the natural frequency of each oscillator and its connection degree in network, i.e., ωi = ki. For model (1), it has been shown that ES could be observed in general network topologies, such as fully connected networks, Erdös-Rényi (ER) networks and scale-free networks of the uncorrelated configuration model19. In the present study, we focus on the case of fully connected networks, which usually are relatively easier to deal with in theoretical treatment when mean-field methods are employed.

Results

The forward phase transition

Typically, an order-parameter is defined as the ensemble average of the complex amplitudes of the oscillators:

where r is the average amplitude characterizing the coherence of the system and ϕ the average phase. Model (1) exhibits first-order synchronization transition with typical FDs, such as triangle, Lorentzian, Gaussian, bimodal Lorentzian, etc. Some examples are illustrated in Fig. 1 and Fig. 2(a).

Figure 1
figure 1

Examples of ES in Eq. (1) with triangle FD (a) and bimodal Lorentzian FD with Δ = 1 and ω0 = 1 (b).

The forward (red lines) and backward (blue lines) transitions in the plots of the order parameter r vs. the coupling strength κ allow defining a clear hysteresis region. The dashed lines are the unstable branch in the backward transitions (see Eq. (31) in the main text). Table I reports the formulae for the various FDs used. In the simulations N = 20, 000 oscillators have been considered and a fourth order Runge-Kutta integration method with time step 0.01 is used.

Full size image
Figure 2
figure 2

(a) Bifurcation diagram and hysteresis area for the first-order PT in Eq. (1) with Lorentzian FD. Besides the incoherent state with r = 0, there are two coherent solutions: the solid blue line r1(κ) is the stable backward PT line, while the dashed blue line r2(κ) is the unstable one. For Lorentzian FD, they are both independent of Δ, i.e., the width of distribution. Red stars, blue circles and green squares are simulation results corresponding to N = 20, 000 and Δ = 1.0, 2.0, 3.0, respectively. (b) Schematic plot of the backward PT in the complex plane. When κ ~ ∞, the stable (solid blue dot) and unstable coherent (solid red dot) solutions are located near the real and imaginary axes, respectively. With the decrease of κ, they move oppositely along the unit circle. Eventually, they collide and disappear at θ = ±π/4, corresponding to κb = 2, in (a).

Full size image

We can analytically solve the critical coupling strength κf for the forward PT in Eq. (1) by using a method similar to that of Refs. 2, 3. The main idea is to consider the continuum limit of the model, i.e., N → ∞. In this case, a density function ρ(θ, ω, t) can be defined, which denotes the fraction of oscillators with frequency ω whose phases have values between θ and θ + dθ at time t. It should be pointed out that ρ(θ, ω, t) is by definition a conditional probability density and cannot be understood as a joint probability density. ρ(θ, ω, t) satisfies the following normalization condition

for all ω and all t. The evolution of ρ(θ, ω, t) is governed by the continuity equation.

There are two ways to solve for the critical points. One is by adding noise to Eq. (1), analyzing the corresponding Fokker-Planck equation and letting the noise strength eventually approach zero. The other is to directly linearize the continuity equation. Remarkably, we obtain the same results from both methods. Therefore, in the following, we will concentrate on the second method, while we refer the interested reader to the Supplementary Information for the details of the Fokker-Planck treatment.

In the continuum limit, Eq. (1) can be written in its mean-field form. Multiplying a factor to both sides of Eq. (2), one has:

Extract the imaginary part from the above equation, i.e., and substitute it into Eq. (1), one gets the dynamical equations in the mean-field form:

where r and ϕ are order parameters defined in Eq. (2). The continuity equation for ρ is:

where the velocity υ is given by

When the coupling strength is relatively small, all oscillators are rotating in the unit circle almost exclusively according to their natural frequencies. This corresponds to the incoherent state ρ0(θ, ω, t) = 1/(2π). One can therefore extract the point at which such an incoherent state loses its stability and thus determine the critical point for forward PT. Specifically, suppose there is a small perturbation from the incoherent state, i.e.,

with . Following the normalization condition Eq. (3), one has that

and the order parameter in continuum form becomes

where

Substitute Eqs. (7), (8) and (11) into Eq. (6), one can obtain its linearized form as

Now η(θ, ω, t) can be further expanded into Fourier series as

where representing the higher Fourier harmonics. By combining Eq. (11) and Eq. (13), one has

i.e., only the first-order harmonic contributes to the order parameter because except n = 0. Similarly, one can have

Adding these two equations, one gets

Substituting Eqs. (13) and (16) into Eq. (12) and comparing the coefficients of term eiθ, one finally obtains

Now, one may seek a type of solution of the form c(ω, t) = b(ω)eλt, where λ is independent of ω. Then the above equation becomes

Equation (18) can be solved in a self-consistent way. Precisely, let . Then b(ω) can be solved from Eq. (18) as: b(ω) = |ω|A/(λ + iω). By substituting this back into the expression of A, one ultimately obtains

Notice that Eq. (19) relates explicitly the coupling strength κ with the eigenvalue λ. If g(ω) is an even function, such as the typical Gaussian and Lorentzian FDs, Eq. (19) can be transformed into

The real part of λ determines the stability of the incoherent state, i.e., when Re[λ] changes from negative to positive, the incoherent state loses its stability. One can then use this condition to determine the critical coupling strength κf for forward PT. For the original Kuramoto model, it has been proven that the eigenvalues are real21. So it is convenient to get κf by directly setting λ = 0 in the corresponding equation of Eq. (20). For the current model, however, there is no guarantee that λ in Eq. (20) is necessarily real. One has to consider the general, thus more difficult, situation where λ is complex rather than real. In fact, it is incorrect to get κf by directly setting λ = 0 in Eq. (20). To solve for λ, one rather has to implement different strategies depending on the specific forms of FD, because in most cases Eq. (20) is a transcendental equation.

Let us now describe the way of overcoming the above difficulty. We take the Lorentzian FD as an example to illustrate the method. By substituting Lorentzian FD (see the formula in Table I) into Eq. (20), one obtains

where x = ω/Δ and z = λ/Δ. The integration of the above equation gives the critical equation relating κ with z, i.e., λ, as

As z is complex and can be written as z = Reiψ, Eq. (22) becomes

Because the RHS of Eq. (23) must be real, then the imaginary part must cancel, i.e.,

At the bifurcation point κf, z (i.e., λ) passes through the imaginary axis and one has cos ψ = 0 and sin ψ = ±1. Using these conditions in Eq. (24), one eventually obtains ln R = 0, i.e., R = 1. In fact, by observing Eq. (24) we find that R = 1 is always its solution regardless of ψ. Now, by substituting z = eiψ back into Eq. (22), one gets

As shown in Fig. 3(a), at the bifurcation point one has , where the subscript denotes bifurcation. Substituting ψb into Eq. (25), we obtain the critical coupling strength for the forward PT as

It should be pointed out that due to the multi-values of complex number, ψb could also take values 2mπ ± π/2 (m = 1, 2, …) at the bifurcation point if we do not limit . This would lead to an infinite number of solutions for transition point. However, those solutions are extraneous mathematically and should be disregarded.

Table 1 Summary of critical equations and critical points. The Table reports the frequency distribution, the critical equation and the obtained value of κf and κb for various FDs: a triangle distribution (first row), a Lorentzian distribution (second row), a Gaussian distribution (third row) and a bimodal Lorentzian distribution (fourth row)
Full size table
Figure 3
figure 3

Schematic illustration of the process through which the incoherent state loses stability.

(a) For Lorentzian FD, the eigenvalue z is on the unit circle in the complex plane. When it crosses the imaginary axis, the incoherent state loses its stability. (b) The solution of the critical coupling strength for the forward PT. From κ = sinc(ψ), we determine κf = sinc(±π/2) = 4.

Full size image

As a bonus of the above solving procedure, one can also determine the stability properties of the incoherent state. As shown in Fig. 3(b), within the interval [0, π], the sinc function monotonically decreases. Therefore, when κ < κf, one finds ψ > ψb = π/2, i.e. z is located at the left of the imaginary axis, indicating that the incoherent state is linearly stable. When instead κ > κf, ψ < ψb = π/2, i.e. z is located at the right of the imaginary axis, implying that the incoherent state is linearly unstable. This argument is also valid for the interval (−π, 0]. Based on these arguments one can conclude that the incoherent state is linearly stable (unstable) before (after) the forward PT critical point in Eq. (1).

The same overall strategy can be actually adopted also for retrieving the critical equation relating κ and λ in the case of other FDs. The steps to be followed are: i) one initially integrates Eq. (20) with specific FDs to obtain the critical equations, as shown in Table I; ii) supposes that λ is complex. As Re[λ] = 0 when κ = κf, one has λ = iy, ; iii) one substitutes λ = iy into the critical equation. Since the RHS must be real, one can then solve for y; iv) one substitutes the proper y back into the critical equation and determines κf. Remarkably, the method turns out to be applicable in many cases and we could solve, using it, for the forward transition points for various typical FDs, including triangle, Gaussian and bimodal Lorentzian. Table I, for instance, is an account of the different FDs used (first column), of the corresponding critical equations (second column, the analogous to Eq. (22) obtained) and of the exact values that can be extracted for κf (third column). We emphasize that, although the current method is successful in dealing with the above typical FDs, it might fail for certain FDs which lead to divergence in the integral of Eq. (20).

The backward phase transition

We now move to solve for the backward PT, starting from the fully coherent state that corresponds to r ~ 1 at κ ~ ∞. For this purpose, the self-consistency method to obtain the critical point κb is implemented1. Suppose the system is in the coherent state for which rκ ≥ 1. It is straightforward to show that, from the mean-field equation, i.e., Eq. (5), the stationary state of the system is

Here, without loss of generality, we have chosen a special coordinate to let the average phase ϕ = 0. Physically, this state consists of two clusters of oscillators with locked phase θp > 0 and θn = −θp < 0, respectively. Now, as long as g(ω) is centered at 0 and symmetric, the oscillators have to split into the two clusters evenly. Thus the order parameter can be greatly simplified as

Using the solutions in Eq. (27), one has

Solving this algebraic equation, one can obtain the two branches of r as

These two coherent solutions only exist when κ ≥ 2. As schematically illustrated in Fig. 2(a), with the decrease of κ the two solutions gradually approach each other and eventually collide with each other and disappear at κ = 2. Therefore, the critical point for the backward PT is κb = 2.

To retrieve the stability properties of the solutions r1 and r2, one notices that in the coherent state all oscillators condensate into two clusters with locked phases shown in Eq. (27). This implies that the original high-dimensional system, i.e., Eq. (1), has drastically collapsed into a two-dimensional one, which is equivalent to a two-oscillator system governed by the following equations

where ω+ > 0 and ω− < 0 are the average frequencies corresponding to those oscillators that are rotating clockwise and counterclockwise, respectively. Then, due to the symmetry of g(ω), one sets ω+ = −ω− = ω and lets Θ = θ+ − θ−, the above equations can be rewritten as

Such an equation has two fixed points when Θ ∈ [0, π], namely arcsin(2/κ) and π − arcsin(2/κ). Through linear stability analysis, one further finds that the former is stable, while the latter is unstable, corresponding to the branches r1 and r2, respectively.

The above analysis provides the following physical picture for the backward PT. As illustrated in Fig. 2(b), when the coupling strength is large enough, i.e. κ ~ ∞, the oscillators evenly split into two clusters that are symmetric with respect to the real axis. There are therefore two solutions whose order parameters correspond to r1 (stable) and r2 (unstable), respectively. The former is very close to the real axis (r1(∞) ~ 1), while the latter is located very close to the imaginary axis (r2(∞) ~ 0). With the decrease of κ, r1 decreases while r2 increases, with the effect of the two solutions gradually approaching the position θp = π/4 and θn = −π/4 along the unit circle and eventually colliding at the critical point κb = 2. The same analysis also reveals the first-order nature of the PT occurring in Eq. (1). Following the mean-field equation, i.e., Eq. (5), the stationary state Eq. (27) is independent of the natural frequency. This implies that with the increase of the coupling strength, when rκ becomes large enough, all oscillators, regardless of their natural frequencies, must simultaneously lock together. In short, the oscillators can only be phase-locked or phase-unlocked, while any partially phase-locked state is strictly forbidden in this system. Therefore, the transition between the incoherent state and the coherent state has no choice but that of being of the first-order type. It has to be remarked that the closed forms of the two solutions (as well as their stability properties) in the backward transition do not depend on the specific FD, but they constitute a generic feature of Eq. (1). Physically, this implies that the dynamics in the coherent state “forgot” its initial configuration, such as initial condition and FD. All our predictions are greatly verified by large scale simulations (N = 20, 000), as shown in Fig. 1, Fig. 2(a) and Table I.

Discussion

In conclusion, the first-order synchronization transition in a generalized Kuramoto model has been analytically studied. We have reported the exact solutions for the critical coupling strengths for both the forward and backward transitions and we carried out the linear stability analysis for both the incoherent and coherent states. In particular, we have demonstrated that the closed forms for both the stable and unstable coherent states for the backward transition are independent on the specific form of FD. The set of results here reported provides insights and a better understanding on the mechanisms leading to the first-order synchronization transition of networked oscillators.