Competing synchronization on random networks

TL;DR

This paper investigates how competing synchronization dynamics in Kuramoto models behave on Erdős–Rényi networks, revealing that network heterogeneity significantly alters phase transitions and states compared to fully connected networks.

Contribution

It extends the analysis of competing Kuramoto models to ER networks, showing the impact of network degree heterogeneity on synchronization phases and transitions.

Findings

High mean degree ER networks exhibit phase diagrams similar to complete networks.

Low mean degree ER networks lack traveling wave states, altering phase diagrams.

Mean-field approximation aligns with simulations only for high mean degree networks.

Abstract

The synchronization pattern of a fully connected competing Kuramoto model with a uniform intrinsic frequency distribution $g(\omega)$ was recently considered. This competing Kuramoto model assigns two coupling constants with opposite signs, $K_1 < 0$ and $K_2 > 0$, to the $1-p$ and $p$ fractions of nodes, respectively. This model has a rich phase diagram that includes incoherent, $\pi$, and traveling wave (TW) phases and a hybrid phase transition with abnormal properties that occurs through an intermediate metastable $\pi$ state. Here, we consider the competing Kuramoto model on Erd\H{o}s--R\'enyi (ER) random networks. Numerical simulations and the mean-field solution based on the annealed network approximation reveal that in this case, when the mean degree of the random networks is large, the features of the phase diagram and transition types are consistent overall with those on completely connected networks. However, when the mean degree is small, the mean-field solution is not consistent with the numerical simulation results; specifically, the TW state does not occur, and thus the phase diagram is changed, owing to the strong heterogeneity of the local environment. By contrast, for the original Kuramoto oscillators, the annealed mean-field solution is consistent with the numerical simulation result for ER networks.

Figures

Equations

$${\dot{\theta }}_i={\omega }_i+\frac{K}{N}\sum _{j=1}^N\sin ({\theta }_j-{\theta }_i),$$

$$\dot{{\theta }_i}={\omega }_i+\frac{K_i}{N}\sum _{j=1}^N\sin ({\theta }_j-{\theta }_i),i=1,2,\cdots ,N,$$

$$Z(t)\equiv R(t)e^{i\psi (t)}\equiv \frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j(t)},$$

$$\dot{{\theta }_i}={\omega }_i+K_{i}R(t)\sin \left(\psi (t)-{\theta }_i\right),i=1,2,\cdots ,N.$$

$$v_1$$

$$\text{vspace}5mm{v}_2$$

$$\text{vspace}5mmv$$

$${\psi }_{\alpha }=\arctan \left[\frac{⟨\sin {\theta }_i{⟩}_{\alpha }}{⟨\cos {\theta }_i{⟩}_{\alpha }}\right]=\arctan \left[\frac{⟨{\theta }_i{⟩}_{\alpha }-⟨{\theta }_i^3{⟩}_{\alpha }/3!+\cdots }{1-⟨{\theta }_i^2{⟩}_{\alpha }/2!+\cdots }\right].$$

$$\dot{\delta }\approx -⟨{\dot{\theta }}_{12}⟩=-(v_1-v_2)=-(p-p_u)(|K_1|+K_2)r_1{r}_2\sin \delta ,$$

$$p_u=\frac{1}{Q+1}.$$

$${\dot{\varphi }}_j={\omega }_j-\Omega -K_{j}R\sin {\varphi }_j.$$

$$R=\frac{1}{N}\sum _j\overline{e^{i{\varphi }_j(t)}},$$

$${\varphi }_j^{\ast }=\arcsin \left[\frac{{\omega }_j-\Omega }{K_{j}R}\right],$$

$$T_j={\int }_0^{2\pi }\frac{d{\varphi }_j}{|{\dot{\varphi }}_j|}=\frac{2\pi }{\sqrt{({\omega }_j-\Omega {)}^2-(K_{j}R{)}^2}}.$$

$$\overline{e^{\mathrm{i}{\varphi }_j}}$$

$$R$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial \theta }\left[\rho (K,\omega ,\theta ,t)\left(\omega +K\frac{Z{e}^{-i\theta }-Z^{\ast }{e}^{i\theta }}{2i}\right)\right]=0.$$

$$\rho =\frac{g(K,\omega )}{2\pi }\left[1+\sum _{n=1}^{\infty }\left(a_{\omega ,K}^n(t)e^{in\theta }+c.c.\right)\right],$$

$${\dot{a}}_{\omega ,K}(t)=-i\omega a_{\omega ,K}+\frac{K}{2}(Z^{\ast }-Z{a}_{\omega ,K}^2).$$

$$Z(t)$$

$$g(K,\omega )=\frac{1}{2\gamma }\Theta (\gamma -|\omega |)\left[(1-p)\delta (K-K_1)+p\delta (K-K_2)\right].$$

$$R=-(1-p)\int_{-|K_{1}|R}^{|K_{1}|R}d\omega g(\omega)\sqrt{1-\Big{(}\frac{\omega}{|K_{1}|R}\Big{)}^{2}}+p\int_{-K_{2}R}^{K_{2}R}d\omega g(\omega)\sqrt{1-\Big{(}\frac{\omega}{K_{2}R}\Big{)}^{2}},$$

$$R=\begin{cases}-\frac{1-p}{2\gamma}\frac{\pi}{2}|K_{1}|R+\frac{p}{2\gamma}\frac{\pi}{2}K_{2}R&{\rm for~{}~{}}R<\frac{\gamma}{|K_{1}|}<\frac{\gamma}{K_{2}}\\[10.0pt] -(1-p)\Big{(}\frac{|K_{1}|R}{2\gamma}\arcsin\frac{\gamma}{|K_{1}|R}+\frac{1}{2}\sqrt{1-\big{(}\frac{\omega}{|K_{1}|r}\big{)}^{2}}\Big{)}+\frac{p}{2\gamma}\frac{\pi}{2}K_{2}R&{\rm for~{}~{}}\frac{\gamma}{|K_{1}|}<R<\frac{\gamma}{K_{2}}\\[10.0pt] -(1-p)\Big{(}\frac{|K_{1}|R}{2\gamma}\arcsin\frac{\gamma}{|K_{1}|R}+\frac{1}{2}\sqrt{1-\big{(}\frac{\omega}{|K_{1}|R}\big{)}^{2}}\Big{)}\\ \hskip 90.00014pt+p\Big{(}\frac{K_{2}R}{2\gamma}\arcsin\frac{\gamma}{K_{2}R}+\frac{1}{2}\sqrt{1-\big{(}\frac{\omega}{K_{2}R}\big{)}^{2}}\Big{)}&{\rm for~{}~{}}\frac{\gamma}{|K_{1}|}<\frac{\gamma}{K_{2}}<R.\end{cases}$$

$$p=\{\begin{array}{ll}\frac{|K_1|+4\gamma /\pi }{|K_1|+K_2}& \mathrm{for}0<R<\frac{\gamma }{K_2},\\ \frac{(\pi |K_1|+4\gamma )R}{\pi |K_1|R+2\gamma \sqrt{1-(\frac{\gamma }{K_{2}R}{)}^2}+2K_{2}R\arcsin (\frac{\gamma }{K_{2}R})}& \mathrm{for}\frac{\gamma }{K_2}<R<\frac{\gamma }{|K_1|},\\ \frac{\gamma \sqrt{1-(\gamma /|K_1|R{)}^2}+|K_1|R\arcsin \left(\frac{\gamma }{|K_1|R}\right)+2\gamma R}{\gamma \sqrt{1-(\frac{\gamma }{|K_1|R}{)}^2}+|K_1|R\arcsin \left(\frac{\gamma }{|K_1|r}\right)+\gamma \sqrt{1-(\frac{\gamma }{K_{2}R}{)}^2}+K_{2}R\arcsin \left(\frac{\gamma }{K_{2}R}\right)}& \mathrm{for}\frac{\gamma }{K_2}<\frac{\gamma }{|K_1|}<R.\end{array}$$

$$p_c$$

$$R_c$$

$$p$$

$$⟨K{⟩}_c=(1-p_c)K_1+p_c{K}_2=\frac{4\gamma }{\pi },$$

$$R=\begin{cases}-\frac{1-p}{2\gamma}\frac{\pi}{2}|K_{1}|R+\frac{p}{2\gamma}\frac{\pi}{2}K_{2}R&{\rm for~{}~{}}R<\frac{\gamma}{|K_{1}|}<\frac{\gamma}{K_{2}}\\[10.0pt] -(1-p)\Big{(}\frac{|K_{1}|R}{2\gamma}\arcsin\frac{\gamma}{|K_{1}|R}+\frac{1}{2}\sqrt{1-\big{(}\frac{\omega}{|K_{1}|r}\big{)}^{2}}\Big{)}+\frac{p}{2\gamma}\frac{\pi}{2}K_{2}R&{\rm for~{}~{}}\frac{\gamma}{|K_{1}|}<R<\frac{\gamma}{K_{2}}\\[10.0pt] -(1-p)\Big{(}\frac{|K_{1}|R}{2\gamma}\arcsin\frac{\gamma}{|K_{1}|R}+\frac{1}{2}\sqrt{1-\big{(}\frac{\omega}{|K_{1}|R}\big{)}^{2}}\Big{)}\\ \hskip 80.00012pt+p\Big{(}\frac{K_{2}R}{2\gamma}\arcsin\frac{\gamma}{K_{2}R}+\frac{1}{2}\sqrt{1-\big{(}\frac{\omega}{K_{2}R}\big{)}^{2}}\Big{)}&{\rm for~{}~{}}\frac{\gamma}{|K_{1}|}<\frac{\gamma}{K_{2}}<R.\end{cases}$$

$$p=\{\begin{array}{ll}\frac{|K_1|+4\gamma /\pi }{|K_1|+K_2}& \mathrm{for}0<R<\frac{\gamma }{|K_1|}\\ \frac{2R\gamma +\gamma \sqrt{1-(\frac{\gamma }{|K_1|R}{)}^2}+|K_1|R\arcsin \left(\frac{\gamma }{|K_1|R}\right)}{\frac{\pi K_{2}R}{2}+\gamma \sqrt{1-(\frac{\gamma }{|K_1|R}{)}^2}+|K_1|R\arcsin \left(\frac{\gamma }{|K_1|R}\right)}& \mathrm{for}\frac{\gamma }{|K_1|}<R<\frac{\gamma }{K_2}\\ \frac{\gamma \sqrt{1-(\frac{\gamma }{|K_1|R}{)}^2}+|K_1|R\arcsin \left(\frac{\gamma }{|K_1|R}\right)+2\gamma R}{\gamma \sqrt{1-(\frac{\gamma }{|K_1|R}{)}^2}+|K_1|R\arcsin \left(\frac{\gamma }{|K_1|R}\right)+\gamma \sqrt{1-(\frac{\gamma }{K_{2}R}{)}^2}+K_{2}R\arcsin \left(\frac{\gamma }{K_{2}R}\right)}& \mathrm{for}\frac{\gamma }{|K_1|}<\frac{\gamma }{K_2}<R.\end{array}$$

$$p$$

$$p_{c,f}=\frac{|K_1|+4\gamma /\pi }{|K_1|+K_2}=\frac{Q+\frac{4\gamma }{\pi K_2}}{Q+1}=p_c,$$

Full text

Keywords: competing Kuramoto model, synchronization, random network

1 Introduction

The Kuramoto model (KM) [1] is a coupled phase oscillator model that exhibits a synchronization transition from an incoherent (IC) to a synchronized state. This model is written as

[TABLE]

where $\theta_{i}$ denotes the phase of oscillator $i$; $\omega_{i}$ is the natural frequency of oscillator $i$, which follows the distribution $g(\omega)$; $K$ is a positive coupling constant; and $N$ is the number of oscillators in the system. We note that each oscillator interacts with all the other oscillators, which is characteristic of a conventional mean-field model. This KM has been used as a prototypical model to understand synchronization phenomena in diverse real-world systems ranging from biological systems such as fireflies flashing to physical systems such as electric power grids [2, 3, 4]. In the original KM, $g(\omega)$ has a Gaussian distribution, and a continuous synchronization transition occurs as the coupling strength $K$ is increased.

The coupling constant is not necessarily positive but can be generalized to a mixture of positive and negative constants when the system contains competing interactions. Here, we consider a KM with competing coupling constants, which is inspired by neural networks with excitatory and inhibitory couplings between neurons, and by spin glass systems with ferromagnetic and antiferromagnetic interactions. For instance, such competing interactions in spin glass systems often result in the frustration of magnetic ordering and the emergence of a glass state. Metastable slow dynamics proceeds across the free energy landscape. Therefore, similar metastable behaviors may be expected to arise in the competing KM ($c$-KM).

The KM with competing interaction was first considered by Daido [5] in 1992. He used the Sherrington–Kirkpatrick random coupling constant $K_{ij}$ between all pairs of oscillators $(i,j)$, which follows a Gaussian distribution with zero mean and a finite standard deviation. Under this type of competing interaction, several interesting behaviors appeared, for instance, quasi-entrainment, diffusive dynamics in phase space, and a slow power-law decay of coherence. Furthermore, the possible existence of a glassy oscillating state, which leads to the so-called volcano transition, was suggested [5, 6].

Further, Hong and Strogatz [7, 8, 9] considered a tractable KM with node-based random couplings $\{K_{i}\}$, in which a $1-p$ fraction of the nodes have a negative coupling constant $K_{1}<0$, and the remaining nodes have a positive coupling constant $K_{2}$. The authors adopted a Lorentzian $g(\omega)$ and used the Ott–Antonsen (OA) ansatz [10, 11]; they found that the synchronized state of this $c$-KM can be characterized by the dynamics of two groups of oscillators roughly separated by an angle $\pi$ in phase space. They are either static or traveling.

The OA method has been successful in the bifurcation analysis of several KMs with unimodal $g(\omega)$ such as a Lorentzian distribution [8, 9, 10, 11, 12, 13]. The reason is that an $N\rightarrow\infty$ continuum of the oscillator degrees of freedom is exactly reduced to a few coupled modes. The Watanabe–Strogatz method is also useful for reducing the degrees of freedom of the coupled nonlinear equations of $N$ identical oscillators [7, 14, 15, 16]. The order parameter dynamics is exactly calculable in these OA and Watanabe–Strogatz reducible systems. However, for a uniform distribution, reduction is hardly achieved. Instead, one can at best seek the stationary solution using the complex self-consistency (SC) equation.

In our previous study [17], we considered a $c$-KM model with a uniform $g(\omega)$ in the interval $[-\gamma,\gamma]$, which is called the competing $c$-Winfree–Pazó ($c$-WP) model. Using the complex SC equation, we showed briefly that a discontinuous synchronization transition with critical behavior, which is often called a hybrid phase transition, occurs. We focused mainly on the two-step discontinuous transitions that arise above the hybrid transition point: the transitions from the IC to the $\pi$ state, and then to the traveling wave (TW) state, where the intermediate $\pi$ state is characterized by a long period of metastability.

In this study, we first recapitulate a general formalism of the complex SC equation that is applicable to general types of $g(\omega,K)$ [17, 18]. Next, we derive explicitly the complex SC equation for the $c$-WP model. In this case, the $\pi$ and TW states can emerge in the synchronized regime. The imaginary part of the SC equation determines the TW solution. A rich phase diagram is obtained, and critical behaviors are derived explicitly as a function of the parameter set $\{p,Q\equiv|K_{1}|/K_{2},\gamma\}$. However, the complex SC equation does not determine the stability of the solution rigorously. Instead, the empirical linear stability of each complex SC equation solution can be analyzed [18]. However, the linear stability did not distinguish the metastable states that arise in the $c$-WP model [17].

Next, we consider the $c$-KM on Erdős–Rényi (ER) networks. First, using the mean-field theory for annealed ER networks, we check analytically whether the synchronization transitions of the all-to-all coupling case occur even on ER networks. We find that the synchronization transitions are similar overall. However, the simulation results show significant differences. On ER networks, the mean degree of the ER network, $\langle q\rangle$, also plays an important role in determining the synchronization transition properties. When $\langle q\rangle$ is small, the TW state can disappear, even though the analytic and numerical solutions for annealed ER networks predict its existence. Moreover, the transition type can change.

This paper is organized as follows. In Sec. 2, we introduce the $c$-KM and investigate the origin of the TW state of the oscillators. The mean-field solution of the TW is also presented. In Sec. 3, we derive the complex SC equation using the original Kuramoto method and the OA ansatz. In Sec. 4, we consider the $c$-WP model and derive the SC equation explicitly. Using the SC equation, we derive the order parameters $R$ and $\Omega$ of the $c$-WP model for various $Q=|K_{1}|/K_{2}$ and the half-width $\gamma$ of the uniformly distributed $g(\omega)$. We obtain rich phase diagrams including IC, $\pi$, and TW phases and continuous, discontinuous, and hybrid phase transitions. In Sec. 5, we consider the $c$-KM model on ER networks. Using the mean-field approach, we obtain the mean-field SC solution for annealed ER networks and the corresponding phase diagrams. Both annealed and quenched ER networks are numerically simulated. The order parameter curves are obtained as a function of $p$ for various values of $\gamma$, $Q$, and the mean degree $\langle q\rangle$. We observe some discrepancies between the mean-field solution and the quenched ER network results, and consider possible reasons for the inconsistencies. In Sec. 6, we summarize our results.

2 The competing Kuramoto model

The $c$-KM on completely connected networks is written as

[TABLE]

where the phase $\theta_{i}$ of each oscillator has an intrinsic frequency $\omega_{i}$, and $\omega_{i}$ follows the distribution function $g(\omega)$. In this section, we do not assume a specific form of $g(\omega)$. $K_{i}$ is the coupling constant of oscillator $i$. A fraction $(1-p)$ of the nodes have a coupling constant $K_{1}<0$, and the remaining nodes have $K_{2}>0$. In addition, $p$ is a control parameter. This $c$-KM can be generalized to an arbitrary number of species of coupling constants.

The synchronization transition is characterized by a complex order parameter at time $t$,

[TABLE]

where $R$ is the coherence of the oscillators, and $\psi$ is the average phase. In the long-time limit, the system falls into the stationary state, in which $Z\approx Re^{i\Omega t}$. The synchronization order parameter $R$ distinguishes the coherent (C) and IC phases. The C phase is further divided into the $\pi$ and TW phases by the TW order parameter $\Omega$. Overall, three phases are possible: the IC ($R=0$), $\pi$ ($R\neq 0,\Omega=0$), and TW ($R\neq 0,\Omega\neq 0$) phases.

Inserting the definition of the complex order parameter into (2) results:

[TABLE]

Depending on the sign of the coupling constant $K_{i}$, the stability of an oscillator at a velocity-balancing position is reversed. Thus, the oscillators are tiered into two groups on the phase circle. Each oscillator with $K_{2}>0$ is attracted toward the average phase $\psi$ with a strength proportional to the coherence $R$, and each oscillator with $K_{1}<0$ is drawn toward the antipode $\psi+\pi$, as shown in the schematic illustration in Fig. 1. Therefore, in the synchronized phase, the oscillators are clustered into two groups, where each positive and negative coupling population forms a separate group. The $K_{1}$ group, i.e., the group of oscillators with the coupling constant $K_{1}<0$, generally shows a broader angular distribution than the $K_{2}$ group, because each $K_{1}$ oscillator repels all the other oscillators, whereas each $K_{2}$ oscillator attracts all the other oscillators.

In the $\pi$ phase, the two groups are balanced at a separation angle $\pi$ and occupy a constant position on the phase circle. In the TW phase, a separation of less than $\pi$ is maintained, where the attractive $K_{2}$ group tries to catch up with the $K_{1}$ group, and the repulsive $K_{1}$ group tries to maintain its distance from the $K_{2}$ group. As a result, the two groups co-rotate around the phase circle at a common angular speed $\Omega$.

2.1 The origin of the TW phase

This non-static steady state is a novel property of $c$-KMs [8]. The TW phase is possible in the competing model owing to the asymmetric interactions. In particular, competing interactions may cause a nontrivial net rotation $\Omega$ in the system. Note that $\Omega$ is different from the average natural frequency, which is zero.

The TW phase can emerge only if the attractive force of the $K_{2}$ group is sufficiently larger than the repulsive force of the $K_{1}$ group, i.e., if $Q\equiv|K_{1}|/K_{2}<1$. Otherwise, a separation of less than $\pi$ is unstable. We construct a mean-field theory below as a first attempt to understand the emergence of the TW phase.

2.2 Mean-field theory for the TW phase

In the TW phase, the system rotates collectively at a nonzero speed. The average angular speed [19] of the groups $K_{1}$ and $K_{2}$ and that of the total system are calculated as follows:

[TABLE]

where $\langle e^{i\theta_{j}}\rangle_{1,2}\equiv r_{1,2}e^{i\psi_{1,2}}$, and the separation between the two groups is denoted as $\psi_{1}-\psi_{2}\equiv\Delta$. For $g(\omega)$ with even symmetry, $\langle\omega_{j}\rangle_{1}=\langle\omega_{j}\rangle_{2}=0$. The angle brackets denote averaging over a group of oscillators. In addition, $\langle\sin\theta_{12}\rangle\equiv\frac{1}{N_{1}N_{2}}\sum_{j\in K_{1}}\sum_{k\in K_{2}}\sin(\theta_{j}-\theta_{k})=r_{1}r_{2}\sin(\psi_{1}-\psi_{2})=r_{1}r_{2}\sin\Delta$. The summation counts only for intergroup interactions because the intragroup interactions cancel out pairwise to zero. When $\Delta=\pi$, $\sin\Delta$ is zero, and therefore $v_{1}=v_{2}=v=0$. However, in the TW phase, $\Delta$ is less than $\pi$, and the angular speed $v$ is nonzero. Now for each group $\alpha=1$ and $2$,

[TABLE]

As a mean-field approximation, we let $\Delta\approx\langle\theta_{12}\rangle\equiv\frac{1}{N_{1}N_{2}}\sum_{j\in 1}\sum_{k\in 2}(\theta_{j}-\theta_{k})$ and $\psi_{\alpha}\approx\langle\theta_{i}\rangle_{\alpha}(1+(\langle\theta_{i}^{2}\rangle_{\alpha}-\langle\theta_{i}\rangle_{\alpha}^{2})/2!)\approx\langle\theta_{i}\rangle_{\alpha}$ for each group; i.e., $r_{1}\approx r_{2}\approx 1$. The approximation is valid as long as the distribution angle of each group is small. Let $\Delta\equiv\pi-\delta$. Assuming that the amplitude components are stable, the phase dynamics of $\delta$ is given as

[TABLE]

where

[TABLE]

Note that for $p>p_{u}$, the $\pi$ state ($\delta=0$) is a stable solution. The $\pi$ state becomes unstable at $p=p_{u}$ as $p$ is decreased, suggesting that a new solution with $\Delta\neq\pi$ (the TW phase) may emerge for $p<p_{u}$. Using the SC equation, we find that the TW solution indeed exists in some interval $[p_{\ell},p_{u}]$. However, this lower bound $p_{\ell}$ is not determined from the mean-field calculation.

3 The complex SC equation

Here, we search for steady-state solutions $Z\approx Re^{i\Omega t}$. For this purpose, we apply the SC method as follows. In the rotating frame $\phi_{j}\equiv\theta_{j}-\Omega t$, each oscillator follows

[TABLE]

Depending on the relative magnitudes of the frequency disorder $\omega_{i}-\Omega$ and the coupling strength $|K_{j}R|$, each oscillator is either phase-locked or drifting. The order parameter is written as

[TABLE]

and the time average, which is denoted by the overline, is taken in the steady-state regime for self-consistency. Each oscillator with a natural frequency satisfying $|\omega_{j}-\Omega|\leq|K_{j}|R$ becomes phase-locked at

[TABLE]

which ultimately contributes to the total synchronization order. In addition, each drifting oscillator satisfying $|\omega_{j}-\Omega|\geq|K_{j}|R$ will repeatedly move ahead of or fall behind the locked population. Their rotation period is calculated as

[TABLE]

Therefore, the contribution of the drifting oscillators to the order parameter $R$ is zero. Indeed, individual drifting oscillators contribute only to the order parameter $\Omega$ as

[TABLE]

After some calculations, the complex SC equation is obtained as follows:

[TABLE]

where $g(K,\omega)$ is the distribution of disorder in the $c$-KM, and ${\rm sgn}(x)=1$ and $-1$ for $x>0$ and $x<0$, respectively.

The imaginary part on the right-hand side of the complex SC equation (14) has to vanish, because $R$ on the left-hand side is real-valued. In the KM, $\Omega$ should generally be given as the mean intrinsic frequency of the system; i.e., the TW order is absent. Therefore, only the real part of the equation must be solved. Indeed, when $g$ is symmetric with respect to $\Omega$, one can easily check that the integrand of the imaginary part becomes odd in $\omega$ and goes to zero. To obtain the TW solution, however, the real and imaginary parts of Eq. (14) must be solved simultaneously for $R$ and $\Omega$.

The complex SC equation can also be obtained by applying the OA ansatz to the continuum version of the KM.

[TABLE]

The density function is assumed to have the form

[TABLE]

where $a_{\omega,K}$ satisfies the equation

[TABLE]

The order parameter is then given as

[TABLE]

In the steady state, Eqs. (16) and (17) also lead to the complex SC equation that we derived as Eq. (14) [18]. For the Lorentzian frequency distribution $g(\omega)=\frac{\gamma/\pi}{\omega^{2}+\gamma^{2}}$, the above integration can be performed in the complex domain, and only a single complex mode, $a^{*}_{-i\gamma,K}$ (for each $K$), contributes to the order parameter. This dimensional reduction does not generally occur for other types of frequency distributions; however, Eq. (17) remains valid. Therefore, exact bifurcation analysis has been applied to limited intrinsic frequency distributions, for instance, the Watanabe–Strogatz transform [14, 15] for identical oscillators or the OA method for Lorentzian-type intrinsic frequency distributions [6, 8, 9, 10, 11, 12, 13, 18, 19]. A general theory for the stability analysis is yet unknown.

Note that the SC method does not determine the stability of the solution in general. The stability of each SC solution must be checked numerically. The empirical linear stability criterion proposed in [18] has been used to identify the stability of the self-consistently obtained solutions in several cases with different intrinsic frequency distributions. However, this criterion is incomplete for metastable states with almost neutral stability for a uniform intrinsic frequency distribution [17].

4 The $c$-WP model

Thus far, we have extended the SC framework to the $c$-KM. Here, we study the $c$-KM with a uniform intrinsic frequency distribution in the range $[-\gamma,\gamma]$ and two coupling constants, $K_{1}<0$ and $K_{2}>0$, at a ratio of $1-p$ to $p$. This model is the $c$-WP model. The distribution $g(K,\omega)$ is written as

[TABLE]

In fact, the KM with a uniform frequency distribution, which we called the Winfree-Pazó (WP) model, is known to exhibit a hybrid synchronization transition, which is a discontinuous phase transition [20] accompanying a singularity [21]. We also observe a hybrid phase transition as $p$ is varied in the $c$-WP model with a critical exponent $\beta_{p}=2/3$ at $p_{c}^{+}$, which is analogous to $\beta_{K}=2/3$ as $K$ is varied at $K_{c}^{+}$ in the WP model [21]. This transition occurs when $Q\equiv|K_{1}|/K_{2}<1$ [Fig. 2(a)]. Furthermore, when $Q>1$, there exists an unstable hybrid transition, as shown in Fig. 2(b).

When $\Omega=0$, i.e., for the IC and $\pi$ states, the imaginary part of Eq. (14) vanishes because the integrand becomes odd in $\omega$. The transition from the IC to the $\pi$ state is characterized by solving the real part of the SC equation,

[TABLE]

where $g(\omega)=\frac{1}{2\gamma}\theta(\gamma-|\omega|)$. In the following, analytic solutions of Eq. (19) are obtained for different parameter values. Note that the integration is evaluated differently depending on the relative sizes of $|K_{1}|R$, $K_{2}R$, and $\gamma$. i) For $Q\equiv|K_{1}|/K_{2}<1$, $\gamma$ can fall in one of the three ranges $(0,|K_{1}|R\,]$, $(|K_{1}|R,K_{2}R\,],$, and $(K_{2}R,\infty)$. ii) For $Q>1$, $\gamma$ can fall in one of the three ranges $(0,K_{2}R\,]$, $(K_{2}R,|K_{1}|R\,],$, and $(\,|K_{1}|R,\infty)$.

4.1 The SC solution for IC and $\pi$

4.1.1 Supercritical hybrid phase transition ($Q<1$)

The integration of Eq. (19) when $Q<1$ yields

[TABLE]

Note that the IC ($R=0$) state is a trivial solution of the above SC equation. The remaining nontrivial solutions correspond to $\pi$ states. After inverting the above equations and solving for $p$, we obtain the inverse function of the order parameter curve $p(R)$ as follows:

[TABLE]

Numerical solutions of the order parameter curve $R(p)$ from the SC equations were presented in [17]. Moreover, the stability of the SC solutions was determined by linear stability analysis [18], which will be discussed later.

The order parameter curve $R(p)$ shows a discontinuous jump of size $R_{c}$ at the critical point $p_{c}$.

[TABLE]

Expanding the intermediate branch $\frac{\gamma}{|K_{1}|}<R<\frac{\gamma}{K_{2}}$ in powers of $\epsilon\equiv(R-R_{c})/{R_{c}}$ after the jump gives

[TABLE]

Therefore, the transition from the IC to the $\pi$ state is hybrid, as $(R-R_{c})\sim(p-p_{c})^{\beta_{p}}$, with $\beta_{p}=2/3$ for $Q<1$. Note that at $p_{c}$, the competing system has a mean coupling strength

[TABLE]

which is equivalent to the critical coupling strength $K_{c}=2/\pi g(0)=4\gamma/\pi$ for the IC-to-C transition in the WP model [20, 21]. Note also the appearance of the critical exponent value for the WP model, $\beta_{K}=2/3$. This is rather natural because the WP model corresponds to a particular case of the $c$-WP model with $p=1$. In addition, successive transitions from the IC to the $\pi$ and then to the TW state are observed in some parameter regimes. Further details are provided in [17].

4.1.2 Subcritical hybrid phase transition ($Q>1$)

The integration of Eq. (19) when $Q>1$ yields

[TABLE]

Similarly, we obtain

[TABLE]

As $R\rightarrow R_{c}={\gamma}/{|K_{1}|}$, for $\epsilon=(R_{c}-R)/R_{c}$:

[TABLE]

As the model approaches the IC-to-$\pi$ transition for $Q>1$, there is an unstable solution that is separated from the IC phase, as shown in Fig. 3. Here, the unstable solution follows $(R-R_{c})\sim(p_{c}-p)^{2/3}$, which gives the same exponent, $\beta_{p}=2/3$; however, $p_{c}$ is approached from the opposite direction $p\to p_{c}^{-}$. In fact, this hybrid transition for $Q>1$ is hidden in the simulations, because it corresponds to an unstable branch. The discontinuous transition for $Q>1$ shows a hysteresis curve that starts at $p_{c,f}=p_{c}$ in the forward direction and at $p_{c,b}$ in the backward direction, where $p_{c,f}$ is given as

[TABLE]

and $p_{c,b}$ and $R_{c,b}$ are determined numerically using the equation

[TABLE]

4.2 Phase diagram

Phase diagrams of the synchronization transitions in the $(p,\gamma)$ plane for various $Q$ values are presented in Fig. 4. In Figs. 4(a) and (b), we consider the case $Q\geq 1$. The phase diagrams contain IC and $\pi$ phases and the hysteresis zone H of the two phases. Both types of dashed lines represent discontinuous transitions; however, the forward transition from the IC to the $\pi$ phase is hybrid. At $Q=1$ in Fig. 4(b), the hysteresis vanishes. The $p=1$ line corresponds to the phase diagram of the WP model, and a hybrid synchronization transition occurs at $\gamma_{h}\approx 0.78$ ($\bullet$) in the WP model [20, 21]. In Figs. 4(c)$-$(e), we consider the case $Q<1$. Part of the $\pi$ state is occupied by the TW phase. The hysteresis region reappears between the IC and TW phases. As the ratio $Q$ is decreased, the interval $[p_{\ell},p_{u}]$ of the TW state becomes broader.

4.3 The empirical stability of the SC solutions

To characterize the stability of the SC solutions, Iatsenko et al. [18] assumed a linear relaxation dynamics of the complex order parameter for small perturbations near each SC solution. The (empirical) stability matrix $\hat{S}$ is reproduced as follows:

[TABLE]

where $F_{R}$ and $F_{\Omega}$ are the real and imaginary parts of the SC complex order parameter, respectively. They are written explicitly as

[TABLE]

In the stability matrix $\hat{S}$, the relative strength of the proportional coefficients of $\delta R$ and $\delta\Omega$ has been determined empirically [18].

The empirical stability of the SC solutions of the $c$-WP model is represented by the blue solid (stable) and dashed (unstable) curves in Fig. 3 and Fig. 2 of [17]. Our numerical result suggests that for neutrally stable solutions, the linear stability criterion of [18] may be partially fulfilled. In our previous study [17], we found neutrally (or weakly) stable $\pi$ solutions, which are stable in the $R$ direction but (nearly) neutral in the $\Omega$ direction. Thus, the order parameter trajectories exhibited some drifting motion near these marginally stable solutions. These solutions characterize metastable states, which capture the system for a long time; however, the system eventually escapes toward the final stable TW solutions.

5 $c$-WP model on ER networks

5.1 The SC equation for ER networks

Now we generalize the $c$-KM to random networks. The Kuramoto equation is written as

[TABLE]

where the normalization $\langle q\rangle$ is equivalent to the expected number of terms in the summation, and $a_{ij}$ is the adjacency matrix. On the ER network, $a_{ij}=1$ if two nodes $i$ and $j$ are connected; otherwise, $a_{ij}=0$. For all-to-all connections, this equation reduces to Eq. (2). Now we consider an annealed ER network, in which

[TABLE]

Here, $q_{i}$ and $q_{j}$ are the degrees of nodes $i$ and $j$, respectively. Each degree follows a Poisson distribution,

[TABLE]

where again $\langle q\rangle$ is the mean degree. When the order parameter is defined as the degree-weighted phase coherence,

[TABLE]

the model is written as

[TABLE]

where we take $\phi=0$ for simplicity. In the stationary limit, the oscillators satisfying $|\omega_{i}-\Omega|\leq|K_{i}|q_{i}H/\langle q\rangle$ are phase-locked, where $\Omega$ denotes the angular velocity of the complex order parameter $He^{i\phi}$. Each locked oscillator contributes a phasor of

[TABLE]

and each drifting oscillator rotating with a period

[TABLE]

makes a time-averaged contribution of

[TABLE]

Thus, the SC equation of the $c$-KM on an ER network is written as

[TABLE]

where we used the shorthand notation $X=\langle q\rangle(\omega-\Omega)/(KH)$.

Eq. (40) can be numerically solved. The corresponding mean-field phase diagram is shown in Fig. 5. The TW regime expands as the network mean degree $\langle q\rangle$ is increased. We confirm that as $\langle q\rangle\rightarrow N$, the solution is reduced to the solution on the all-to-all network shown in Fig. 4(c). The stability of these mean-field solutions was tested by numerical simulations on the annealed ER network governed by Eq. (36). The black symbols in Fig. 6 represent the stationary states of the simulation, which were reached after waiting a sufficiently long time. The solid (dotted) blue curves in Fig. 6 denote stable (unstable) solutions of the SC equations under the linear empirical stability criterion, which was calculated by a method similar to one in the literature [17, 18]. We observe that the TW and $\pi$ phases also appear in the annealed ER network. In addition, the empirical stability condition applies in this mean-field solution in most of the regime; however as before, metastability does arise in the $\pi$ states (solid blue), as they are covered by the simulation data points only for an intermediate range of time. Therefore, the mean-field theory is as successful for annealed ER networks as the SC theory is for all-to-all networks in terms of investigating the synchronization phase transition in the $c$-WP model.

5.2 Synchronization transitions of $c$-WP model on ER networks

Next, we compare the results of numerical simulations of $c$-KM on quenched ER random networks, which use the conventional adjacency matrix (32). In Fig. 6, the simulation data points of networks starting from a C/IC state are denoted by green ($\blacktriangle$)/red ($\bullet$) symbols. A comparison with the mean-field results on the annealed ER networks reveals similarities and differences. The TW and $\pi$ phases also appear in the $c$-KM model on ER networks, but only for a sufficiently large mean degree $\langle q\rangle$. Moreover, the order of transition may depend on $\langle q\rangle$. In the ER networks with low mean degree [Figs. 6(a)$-$(d)], the TW is not observed. The synchronization transition to the $\pi$ phase requires a fraction $p$ of $K_{2}$ oscillators that is larger than $p_{c}\approx 0.66$, and the transition is continuous; however, the $p_{c}$ value in the mean-field result is smaller, and the transition can be discontinuous. For an intermediate mean degree [Figs. 6(e)$-$(h)], we observe the emergence of a TW phase, as in the mean-field prediction; however, the order of transition depends on the frequency half-width $\gamma$. In Fig. 6(e), we unexpectedly find that $p$ exhibits hysteresis behavior in the range $[0.54,0.55]$, which did not occur in the mean-field calculations and simulations. Note that in Fig. 6(g), the TW regime appears, although the angular velocity is very small. For higher mean degree [Figs. 6(i)$-$(l)], the hybrid jump transition at $p_{c}$ is gradually restored.

The mean-field theory of the $c$-WP model works for the annealed ER networks (33). Note, however, that the annealed ER networks are weighted and fully connected. By contrast, the synchronization phase transition of the $c$-WP model on the quenched ER networks differs significantly, especially when the mean degree $\langle q\rangle$ is small. In the $c$-WP model, the annealed network approximation fails for the quenched ER network, in contrast to the case of ordinary KMs on complex networks [22]. The main reason for this failure is the heterogeneity of the local environment, which is caused by the network disorder and the oscillators with a negative coupling constant $K_{1}$. In particular, for a small mean degree $\langle q\rangle$, the number of connections to the neighboring $K_{2}$ nodes is even smaller. Hence, a local order parameter is subjected to large finite-$\langle q\rangle$ fluctuations, and the sign of $K$ is dominant in the competing synchronization.

To investigate this effect, we first measured the link coherence $R_{ij}\equiv a_{ij}\left|\frac{1}{T}\int_{0}^{T}e^{i[\theta_{i}(t)-\theta_{j}(t)]}dt\right|$ [23], as shown in Fig. 7. For a sparse connection, $\langle q\rangle=10$ [Fig. 7(a)], the link coherence matrix can be blocked; for the intraspecies interaction between $K_{1}$ and $K_{1}$ oscillators and $K_{2}$ and $K_{2}$ oscillators, coherences are high, and for the interspecies interaction between $K_{1}$ and $K_{2}$ oscillators, coherence is low. In addition, the values within the block are random and highly heterogeneous. The link coherence of the $K_{1}$ and $K_{2}$ oscillators is distributed within a broad range of values. Such large fluctuations can cause the mean-field theory to break down. As the mean number of neighbors $\langle q\rangle$ increases, the border between $K$-species gradually fades, and the link coherence among $K_{2}$ oscillators eventually becomes homogeneous. However, noisy features remain even when $\langle q\rangle=160$. By contrast, the link coherence of the fully connected, annealed ER network exhibits better organization, as shown in Fig. 7(d). Note that the link coherence has a homogeneous value over an extended area, including all of the $K_{2}$ population and some of the $K_{1}$ population on the high-frequency sides. The corresponding oscillators have been entrained into a TW cluster. In addition, the low-coherence striped patterns correspond to some $K_{1}$ oscillators that are very weakly entrained with the TW cluster. For those oscillators, the density of connections $q_{i}$ is small and/or the natural frequency values are distant from the TW velocity. They remain detrained from the TW cluster.

Finally, we looked at snapshots of the oscillators’ phases and phase velocities [Fig. 8] and the time-averaged mean velocity $V_{i}\equiv\frac{1}{T}\int_{0}^{T}\dot{\theta}_{i}(t)dt$ versus local order parameter $R_{i}\equiv\frac{1}{T}\int_{0}^{T}\left|\frac{1}{q_{i}}\sum_{j}a_{ij}e^{i\theta_{j}(t)}\right|dt$ [Fig. 9]. Note that $\langle q\rangle=10$ has large fluctuations, and it cannot divide the $K_{1}$ and $K_{2}$ oscillators into two groups [Fig. 8(a)]. As the network becomes more connected, the clusters of $K_{1}$ and $K_{2}$ oscillators become much clearer, as shown in Figs. 8(b), (c), and (d). In Fig. 9(a), the TW regime is absent. In Figs. 9(b) and (c), the TW order becomes clearer, and the fluctuations in the local order parameter $R_{i}$ gradually decrease. Note that in the annealed network [Fig. 9(d)], the local order parameter has become homogeneous on all the oscillators. Hence, for the annealed network, the mean-field theory becomes valid. From Figs. 7 and 9, we conclude that as the network becomes more connected, the homogeneity of the local environment is eventually restored.

6 Summary and discussion

We investigated the synchronization transitions in the $c$-KM with a uniform intrinsic frequency distribution on completely connected networks and on ER networks. For the former case, we presented the SC solutions and the critical exponent of the IC-to-$\pi$ hybrid synchronization transition explicitly. Phase diagrams including the IC, $\pi$, and TW phases were obtained analytically for different ratios $Q\equiv|K_{1}|/K_{2}$ of the coupling constants as a function of $p$ and $\gamma$, fractions of nodes with positive coupling constant $K_{2}$, and half-widths of the uniform distribution, respectively. Various synchronization transition types and critical behaviors were identified analytically or numerically. The linear stability of each SC solution was tested by numerical simulations.

Further, the synchronization transition was investigated on the random network. Some of the synchronization features were found to depend on the network mean degree. When the mean degree is large, the behaviors of the synchronization transition are similar overall to those on the completely connected networks. However, when the mean degree is small, the TW phase was not observed in the numerical simulations, even though the mean-field solution on the annealed ER networks predicts the existence of the TW state when the other control parameters are suitable. Transition features of all-to-all networks and ER networks appear differently. We note that the mean-field theory cannot provide a complete understanding of the competing synchronization transition on the ER networks and that advanced mean-field theories are necessary.

This work was supported by the National Research Foundation of Korea by Grant No. NRF-2014R1A3A2069005.