Exact explosive synchronization transitions in Kuramoto oscillators with time-delayed coupling

Synchronization commonly occurs in many natural and man-made systems, from neurons in the brain to cardiac cells to power grids to Josephson junction arrays. Transitions to or out of synchrony for coupled oscillators depend on several factors, such as individual frequencies, coupling, interaction time delays and network structure-function relation. Here, using a generalized Kuramoto model of time-delay coupled phase oscillators with frequency-weighted coupling, we study the stability of incoherent and coherent states and the transitions to or out of explosive (abrupt, first-order like) phase synchronization. We analytically derive the exact formulas for the critical coupling strengths at different time delays in both directions of increasing (forward) and decreasing (backward) coupling strengths. We find that time-delay does not affect the transition for the backward direction but can shift the transition for the forward direction of increasing coupling strength. These results provide valuable insights into our understanding of dynamical mechanisms for explosive synchronization in presence of often unavoidable time delays present in many physical and biological systems.

In many physical, biological and technological oscillatory systems, useful function emerges from collective synchronization of an ensemble of constituent oscillators. Examples include working of neurons in the brain 1-3 , phase-locking of Josephson junction arrays 4,5 , the dynamics of power grids 6 . The Kuramoto model 7,8 , originally formulated to simplify the Winfree's coupled oscillator model for the circadian rhythms of plants and animals 9 , remarkably generalizes to explain phase synchronization phenomena in these examples and many more 10,11 . Transitions to or out of synchronization as a consequence of changing coupling strength are analogous to phase transitions studied in statistical physics such as ferromagnetic, superconductive and thermodynamic transitions 12 .
Collective synchronization of coupled phase oscillators depends on several factors: intrinsic frequency distribution, coupling strength, interaction time-delays, network and structure (coupling strength or topology)dynamics (frequency) relation. As the coupling strength is changed across a certain critical value, the transition from incoherent to coherent, or coherent to incoherent states takes place smoothly (the second-order phase transition like) or abruptly (the first-order like). When coupling is associated with oscillator characteristics or outputs, abrupt transitions to synchrony can occur with hysteresis in a variety of coupled oscillator systems, including in Josephon junction arrays 13 , in complex networks of oscillators 12,[14][15][16][17][18] , and in frequency-weighted, mean-field coupled system of Kuramoto models 19 . Time-delay in mean-field coupling can also make the synchronization transition abrupt 20,21 . Despite our current advanced understanding of synchronization transitions in a variety of these systems, the effects of time-delay and frequency correlated mean-field coupling (structure-dynamics relation) in phase synchronization remain to be explored.
Interaction time-delays and structure-function interdependence are usually unavoidable characteristics of spatially distributed, adaptive oscillatory systems, such as neurons in the brain that has a functional organization 3,22,23 . In such systems, smooth or abrupt synchronization transitions may help us to distinguish between normal and abnormal functioning, such as pconerceptual decision-making 24 as an example of normal functions and epileptic seizure as dysfunction 18,25 .
In this work, we analyze a generalized Kuramoto model of time-delay coupled phase oscillators with frequency-weighted global coupling for stability of incoherent states and coherent states and derive the exact analytical solutions for the critical coupling strengths at different time delays in both directions of increasing (forward) and decreasing (backward) coupling strengths. Here, as a general result, we will come to show that the time delay coupling can affect the abrupt synchronization transition only in the forward direction and not in the backward direction in various network topologies and distributed time-delays.

Methods and Results
We consider the following generalized Kuramoto model with time-delay and frequency-weighted coupling: Here, the coupled system consists of N number of oscillators, each with θ i (t) as the instantaneous phase at time t, θ  t ( ) i its derivative and ω i as the natural frequency. A set of N natural frequencies is drawn from a zero-centered symmetric distribution function (g(ω) = g(−ω)). The coupling strength k is modulated by |ω i |. The heterogeneity of couplings thus achieved can represent characteristics of adaptive oscillator systems commonly found in nature. The heterogeneity of couplings thus achieved can represent characteristics of many functionally organized and spatially distributed oscillator systems, some relevant examples of which include power grid networks 12,26 , social communication networks 27,28 and brain neuronal oscillatory networks 3,29 .
Here, we consider the case of fully connected networks so that we can use the mean-field approach and the continuity equation for the time-evolution of instantaneous phase distribution ρ(θ, ω, t) on a unit circle. We define an order parameter r for t − τ time by where r characterizes phase coherence and φ the average phase of the coupled system. In stationary state, the definition of Eq.
(2) will be equivalent to the traditional definition of The imaginary part of Eq. (3) is: With Eqs (1) and (4), we obtain: Dependent on the coupling strength κ and time delay τ for a given frequency distribution g(ω), the coupled system as represented in Eq. (1), can show phase coherence (r > 0), or incoherence (r ≈ 0). For the forward (incoherent to coherent state) transition, we linearize the continuity equation around the incoherent state (ρ(θ, ω, t) = 1/2π) and obtain the critical coupling strength K f for the forward direction. For the backward (coherent to incoherent state) transition, we start with fully coherent state (r = 1) at sufficiently large κ and use the self-consistency approach on the main mean-field equation to obtain the critical coupling strength K b in the backward direction. Here, we show our calculations for a zero-centered Lorentzian distribution of frequencies, but the calculation method can be applied to any smooth symmetric frequency distribution.

N
, the probability density function (ρ(θ, ω, t)) that represents the fraction of oscillators with frequency ω whose phases are distributed between θ and θ + dθ satisfies (i) the normalizing condition: ∫ ρ θ ω = π t ( , , ) 1 0 2 and (ii) the incoherent state value ρ 0 (θ, ω, t) = 1/2π uniformly distributed over the unit circle. We introduce a small perturbation to a completely incoherent state: The continuity equation for ρ is: The flow velocity function v(t) is By sustituting (6), (8), (10) into (9) and considering the consistency of O(ε) terms, we have Here η(θ, ω, t) can be expanded into the following complex Fourier series: i i η ⊥ represents higher Fourier harmonics terms. Now, (13) and (14), i i Using (12), (15) into (11) and comparing the coefficients of e iθ : We now look for a separable solution: c(ω, t) = b(ω)e λt in the equation (16): This leads to: With for a standard Lorentzian distribution, we have: As R = |λ| > 0, κ passes through the bifurcation point when λ = iR or λ = −iR. When λ = iR, ln (20) When λ = −iR, We set κ κ = 1 2 (complex conjugate pairs) to seek for real critical values. When τ and R satisfies: κ 1 = κ 2 and both are real. Equation (22) has many number of intersection points. Only one R = R 0 > 0 is a unique efficient solution and all others are extra roots. The forward critical value of κ is determined by the value of R 0 as follows: All unique solutions in (23) are greater than or equal to 4 for any τ ≥ 0, as shown in Fig. 1(a).
Backward phase transition. Set a rotating frame with the average phase of the system, Here 〈ω〉 is the average frequency of the oscillators. For an even symmetric distribution of g(ω), 〈ω〉 = 0.
, the mean field equation (5) can be transferred into: In the coherent state, all the oscillators are phase locked.  [14], which helps to determine the critical coupling to be K b = 2 at r(κ) ≈ 0.7 for the forward direction. Δθ i+ and Δθ i− represents the two groups in the equation (26).
i i The above equation (29) has four possible solutions of r(κ) for different κ, all of which are complex for κ < 2, two of which are positive for κ ≥ 2, and only one of which increases for increasing κ. Thus, the equation for a viable solution is: κ = κ κ κ κ , for which the first solution with zero imaginary part occurs at κ = 2. Hence, the backward critical value of κ becomes K b = 2 for any even, symmetric distribution function g(ω) (Fig. 1(b)). The inset in Fig. 1(b) shows the real (blue line) and imaginary (red dashed line) parts of r(κ) at different κ: the imaginary part becomes zero at κ = 2 and remains so for increasing κ and, at κ = 2, r(κ) ≈ 0.7 (predicted value of r at the backward transition point). Fig. 2, we numerically verify the above analytical results in a case of frequency-weighted, fixed time-delayed all-to-all coupling. We show that these findings from all-to-all, fixed time-delayed networks hold also true for the cases of sparsely connected networks (Fig. 3(a-d)) and with distributed time-delayed couplings ( Fig. 4(a-d)). We now turn to extend the above results of fully connected networks to a more realistic case of sparsely connected networks. For that, we consider the random Erdos-Renyi (ER) network as an example. For an uncorrelated network, we follow refs 16,30 to rewrite Eq. (1) as

Numerical results and extensions. As shown in
where P(k), 〈k〉, ρ(k; θ, t) represent the degree distribution, average degree, and density of the nodes with phase θ at time t for a given degree k, respectively, and the term h(t) takes into account time fluctuations and is given by , where "Im" stands for the imaginary part. In the thermodynamic limit, the term h(t) can be neglected when the average degree 〈k〉 is large enough 30 . Consequently, Eq. (2) can be rewritten as Neglecting the fluctuation h(t) and substituting Eq. (31) into Eq. (30), we obtain This is exactly Eq. (5). Therefore, the results obtained from Eq. (5) should also work for the case of non-fully connected networks, provided that the term h(t) in Eq. (30) can be neglected. In numerical simulations, the ER networks with a larger average degree 〈k〉 will satisfy the mean-field approximation (30-32) better. Figure 3(a-d) show the dependence of r on κ for different average degrees 〈k〉 and time delay τ, respectively, with (a) 〈k〉 = 10 and τ/2π = 3.2 × 10 −4 ; (b) 〈k〉 = 100 and τ/2π = 3.2 × 10 −4 ; (c) 〈k〉 = 10 and τ/2π = 0.16; and (d) 〈k〉 = 100 and τ/2π = 0.16. Comparing the case of 〈k〉 = 100 in Fig. 3(b,d) with that of 〈k〉 = 10 in Fig. 3(a,c), we see that the former has a larger loop than the latter. Then, if we compare the case of 〈k〉 = 100 in Fig. 3(b,d) with that of fully connected network (〈k〉 = 999) in Fig. 2(b,c), we find that latter has a larger loop than the former. Thus, the size of loop will monotonously decrease with the decrease of 〈k〉. On the other hand, comparing the case of τ/2π = 3.2 × 10 −4 in Fig. 3(a,b) with that of τ/2π = 0.16 in Fig. 3(c,d), we see that their differences are not significant, indicating that the explosive synchronization is robust to the time delay τ. In sum, Fig. 3(a-d) tells us that for different 〈k〉 and τ, there is always a hysteresis loop and the backward K b is always close to κ = 2, confirming the theoretical extension of Eq. (30).
We now extend the findings to the case of non-uniform delays. For this purpose, we let each node have a different τ i and let each τ i be taken from a random uniform distribution with the average 〈τ〉/2π = 0.16 and standard deviation σ. To focus on the effect of distributed τ i , we take the fully connected network as an example. Figure 4(a-d) show the dependence of r on κ for different σ, respectively, with (a) σ = 0.001; (b) σ = 0.01; (c) σ = 0.1; and (d) σ = 1.0. We find that the backward K b is always close to κ = 2 and the forward K f is only slightly different for different σ, indicating that the explosive synchronization is robust to the distribution of time delay.

Discussion and Conclusions
Here, we have generalized the Kuramoto model to include fully connected, time-delayed, frequency-weighted coupling and analytically derived the exact formulas for critical transitions to or out of explosive (abrupt, first-order like) phase synchronization and extended these results to sparsely connected networks and distributed time-delays. We used |ω|-weighted coupling previously used in a non-delayed system 17 . This scheme is one of the ways to obtain and maintain explosive transitions to or out of synchronization. The frequency-based weighting scheme is relevant to many functionally organized and spatially distributed oscillator systems, some relevant examples of which include power grid networks 12,26 , social communication networks 27,28 and brain neuronal oscillatory networks 3,29 . In the example of power grid network, a network consists of Kuramoto oscillators, where the weighted coupling coefficient between two oscillators is related to their own natural frequencies 26,27 . In the example of communication networks, an extrovert contacts his or her neighbors more frequently than an introvert. If we define the contact between two individuals as a kind of coupling and the frequency of contacts as coupling strength, the coupling strength becomes correlated with the characteristics of individuals, i.e., a kind of natural frequency of human interactions 27 .
In addition to the frequency-weighting scheme, there are other documented ways to induce and maintain explosive synchronization, which are: (i) uniform frequency distribution in all-to-all network topology 7,8,14 , (ii) time-delayed coupling 20 , (iii) order parameter-dependent coupling 13 , (iv) scale-free network topology and correlation between intrinsic frequencies and node degre 15 , and (v) coupling based on weighting procedure with network link frequency mismatch and link betweeness 31 . Taken these cases together, it leads to a general notion that additional constraints or 'inertia' in the system, such as avoiding close frequencies, having time-dependence in parameters or network structure -dynamics correlation, can equally induce or enhance explosive synchronization.
In sum, we generalize the Kuramoto model of globally coupled phase oscillators with time-delay and oscillation frequency-modulated coupling considering its relevance to adaptive physical, biological or technoligcal oscillators. We have analytically and numerically studied the stability of first-order synchronization in this generalized Kuramoto model. We have found the exact formulas for the critical coupling strengths at different time delays in both the increasing (forward) and decreasing (backward) directions of coupling strengths. We find that time-delay does not change the transition in the backward direction but can shift the transition for the forward direction. These results, consistent across sparsely connected networks and networks with distributed time delays, provide useful insights into our understanding of dynamical mechanisms leading to explosive synchronization in presence of often unavoidable time delays in realistic spatially distributed and functionally organized systems. We envision that our theoretical work may encourage future research on abrupt collective synchronization in models