First-order like phase transition induced by quenched coupling disorder

Synchronization is a collective dynamical phenomenon that manifests itself in an impressive range of physical and biological systems.^4,5 Here, we explore the emergent collective behavior exhibited by a population of all-to-all coupled oscillators with random symmetric coupling strengths, either positively or negatively valued, and subject to thermal noise. This system is closely related to the mean-field $XY$ model with random coupling disorder, also known as the spin-glass model;^6,7 however, interactions in the $XY$ model are asymmetric rather than symmetric, and coupling strengths obey a normal distribution, rather than being randomly drawn from two distinct values. Furthermore, in the research fields of coupled oscillator networks and synchronization theory, much attention has been paid to “oscillator glasses”; i.e., the possibility of glass-like behavior has been explored by investigating the relaxation dynamics of such systems.^8–14 In this context, we note that for analytical tractability, most studies on oscillator glasses focused on “separable” disorder in the coupling, where the interaction between two oscillators is given by the product of two coupling strengths associated with the adjacent oscillator nodes individually; e.g., $Kij=Jsisj$⁠, where $si(j)=±1$ and $J$ = const.

In this study, we were instead interested in the behavior induced by non-separable symmetric coupling interactions between oscillator nodes whose interaction strengths were randomly chosen from either of two values, one negative (⁠$Kn<0$⁠) and positive (⁠$Kp>0$⁠). In exploring the collective behavior of the system, we focused on determining the critical threshold for the transition from incoherent to coherent oscillatory behavior, especially in the possibility of a universal character regarding models with varying types of coupling interactions. Specifically, we compared the synchronization threshold of the present model (random non-separable coupling interactions) with the one found for models studied previously by one of the authors, where non-separable interactions are characterized by uniform input or output strengths between oscillator nodes.^15–17 In doing so, we considered the effects of quenched random coupling disorder and of thermal noise whose strength is controlled by a temperature $T$⁠.

This paper is structured as follows. In Sec. II, we introduce the model with its governing equations. Section III discusses the emergent collective behavior of the system for zero noise (⁠$T=0$⁠) and explains the phase transition from the incoherent to coherent oscillator state. We further investigate this transition via linear stability analysis. In Sec. IV, we explore the effect of non-zero thermal noise (⁠$T>0$⁠) on the collective behavior of the system based on numerical simulations and a mathematically rigorous mean-field theory.^1,3 This article concluded with a brief summary and discussion in Sec. V.

We consider $N$ coupled oscillators $i=1,…,N=:[N]$ whose phases $θi(t)∈R/2πZ$ obey the dynamics given by

where $Kij$ is the coupling strength between oscillator nodes $i$ and $j$⁠, drawn from the bimodal distribution function

where positive (⁠$Kp>0$⁠) and negative (⁠$Kn<0$⁠) coupling occur with probability $p$ and $1−p$⁠, respectively. We require that the interactions be symmetric between oscillator nodes, i.e., $Kij=Kji$ for all $i$ and $j$⁠; the model definition implies vanishing self-interaction, which is equivalent to letting $Kii=0$ for all oscillators $i∈[N]$⁠. Rescaling time, $t↦Kpt$⁠, we may reduce the number of effective system parameters by introducing the positive valued coupling ratio$Q:=−Kn/Kp>0$⁠. For now, we shall refer to the unscaled coupling, $Kij$⁠, to complete a physical heuristic discussion until Sec. III A; from Sec. III B onward, we shall make use only of the rescaled version of the coupling to perform various analyses (i.e., let $Kp=1$ and $Kn=−Q$⁠).

Note that it is possible to write out the coupling matrix explicitly in terms of a linear combination of a constant and the Erdösz-Rényi (ER) network, which will prove useful in the subsequent analysis,

where the ER network is defined via the symmetric adjacency matrix ${AijER}i,j∈[N]$ with $AiiER=0,i∈[N]$⁠, taking on the value 1 with probability $p$ and 0 otherwise; i.e., let $AijER=AjiER$ be independent Bernoulli variables, $P(AijER)=p$ with $0<p<1$ for $1≤i<j≤N$⁠, and $AijER=0$ otherwise.

The presence of attractive (⁠$Kij>0$⁠) vs repulsive coupling (⁠$Kij<0$⁠) side-by-side may seem unusual for physical systems, but they are common in the context of excitatory/inhibitory interactions occurring (neuro-) biological systems^18,19 and of social systems²⁰ consisting of agents with conformist and contrarian behavior;^15–17 however, while those studies considered coupling with uniform input strength (⁠$Kij=Ki$ for all $i$⁠) or uniform output strength (⁠$Kij=Kj$ for all $j$⁠) from each oscillator node, here, we are concerned with non-uniform (but symmetric) interactions similar to spin-glass models.^6,7 The last term $ηi(t)$ in Eq. (1) models thermal noise in the system, i.e., mean-centered $(⟨ξi(t)⟩=0)$ Gaussian noise with variance

where $T$ defines the temperature via the Einstein–Smoluchowski relation $D=kBT/γ$⁠, where $kB$ is the Boltzmann constant $kB$ and we set the fraction coefficient to $γ=1$⁠.

The presence of both attractive and repulsive interaction induces frustration in the oscillator system. In this context, we note that Eq. (1) is related to the mean-field $XY$ spin-glass model⁷ with random coupling disorder. Introducing the Hamiltonian

the Langevin dynamics in Eq. (1) can be expressed as $θ˙i=−∂H∂θi+ηi(t)$⁠. This formulation allows the interpretation of the model as the “overdamped” version of the Hamiltonian dynamics at finite temperature $T$⁠; i.e., the inertia term proportional to $θ¨i$ is negligible compared to the damping term proportional to $θ˙i$⁠.

To examine how random coupling affects the collective behavior, we first explore the case of zero temperature (⁠$T=0$⁠) to understand which collective states are possible in the system (1). To gain such insight, consider two extreme cases in $Kij$ in Eq. (2):

• For $p=1$⁠, we have uniform coupling with $Kij=Kp(>0)$ for all $i$ and $j$⁠, which can be interpreted as ferromagnetic coupling. The only energetically favorable state of the system is then fully coherent (ordered) since it possesses the lowest energy with identical phase angles (⁠$θi=θj$ for all $i$ and $j$⁠).

• For $p=0$⁠, the coupling is uniform with antiferromagnetic character, $Kij=Kn(<0)$ for all $i$ and $j$⁠. The system then exhibits the incoherent state with random phases $θi∈[0,2π)$⁠, which corresponds to zero energy in the average sense as $N→∞$⁠.

In between these two cases, we expect a phase transition from the incoherent state to the fully coherent state at a certain critical value $p=pc$⁠. This hypothesis is corroborated by the following two theoretical arguments, which predict the critical threshold, $pc$⁠, for the transition from incoherence to coherence (disorder to order).

The transition point $pc$ can be directly deduced from the properties of the distribution $h$⁠. Given that the random coupling $Kij$ is the only source of disorder in the system, we expect that the critical $pc$ is controlled by its mean value, $⟨Kij⟩=pKp+(1−p)Kn$⁠, for the distribution given by Eq. (2). Specifically, we expect the transition to occur when $⟨Kij⟩$ changes sign, i.e., when the coupling switches from being predominantly negative to being predominantly positive. This is the case when $⟨Kij⟩=0$⁠, which implies that

The plausibility of this result becomes evident when considering the special case $Q=1$⁠, where negative and positive coupling is equally distributed: attractive coupling dominates over repulsive coupling exactly when $pc=1/2$ and thus the onset of coherence is induced.

Alternatively, the threshold $pc$ can be obtained by an energetic argument. For the fully coherent state, the $XY$ oscillators possess the same phases, $θi(t)=θj(t)$ for all $i$ and $j$⁠, which implies that $R=1$ in Eq. (10). Accordingly, the energy density $ε:=H/Ntot$ (energy per spin/oscillator) for the fully coherent state in an average sense can be stated as

where $Np$ (or $Nn$⁠) is the number of positive (or negative) couplings among the total number of couplings $Ntot=N(N−1)/2$⁠; in the thermodynamic limit (⁠$N→∞$⁠), the energy density becomes

It is clear that by definition, $Ecoh≤0$ for all values of $0≤p≤1$ and $Q>0$⁠. Note that for this special case of the fully coherent state, the energy relates directly to the average coupling strength, $Ecoh=−⟨Kij⟩$⁠.

Similarly, the energy density for the incoherent state in the thermodynamic limit can be expressed as

where the (ensemble) averages of phase differences are defined as $ξ:=⟨cos⁡(θi−θj)⟩$ and $ξ′:=⟨cos⁡(θi′−θj′)⟩$⁠. Since all oscillators in the incoherent state assume arbitrary phase values over time with equal probability, we have that $ξ=ξ′=0$⁠, and therefore, $Einc=0$ for $N→∞$⁠. At the transition (⁠$p=pc$⁠), continuity of the energy density demands that the two energy densities match, $Ecoh=Einc$⁠, which implies $(Kp−Kn)pc+Kn=0$⁠. As a consequence, we obtain the same result as in (6).

To see that such a phase transition is indeed possible, and to confirm our theoretical prediction for the critical threshold (6), we numerically integrated Eq. (1) for the time $Mt=2×105$ using the Heun method²¹ with discrete time increments of $Δt=0.01$⁠. Initial phases were drawn from the i.i.d. uniform distribution on ${θi(0)}i∈[N]∈[−π,π)$⁠. Transient behavior possibly occurring during the initial time $t<3/4Mt=1.5×105$ was discarded until an equilibrium state was reached, and observed quantities were averaged over the remaining simulation time. We measured the collective behavior of the system using the complex order parameter

where $|⋅|$ denotes the absolute value and $⟨⋅⟩t$ the average over the time interval $3/4Mt≤t≤Mt$⁠. In the framework of the $XY$ spin model (5), the quantity $R$ corresponds to the magnetization of the material for any given spin configuration.

Figure 1 displays the behavior of $R$ as a function of $p$ with $Q=1$ fixed while varying the system size $N$⁠. We observe a phase transition from the incoherent state near $R=0$ to the fully coherent state near $R=1$⁠. Our simulations indicate a critical transition at $p=pc=1/2$⁠, which agrees with our theoretical arguments further above (6). The slope of the curve $R=R(p)$ increases for larger system sizes near the transition. This suggests the presence of a first-order like transition at the threshold $pc$⁠, as one would expect for the thermodynamic limit $N→∞$ based on the argument made further above [see Eq. (6) and the inset in Fig. 1]. The absence of a perfectly discontinuous transition, characteristic of a first-order transition, is explained by finite size effects and transient dynamics near the transition that smooth the transition.

The behavior of $R$ as a function of $Q$ is studied in Fig. 2. We observed that the critical value $pc$ grows with increasing coupling ratio $Q$⁠. Augmenting the value of $Q$ increases the ratio of the negative coupling strength $Kn$ vs positive coupling strength $Kp$⁠. It is clear that, as a result, a larger value of $p$ is required to reach the fully coherent state. We numerically estimated the critical value for $pc$ based on the simulations using two methods: (i) We measured the value $p=pc$ for which $R(t)>ε$⁠, where we chose $ε=0.02$⁠. [Since the value of $R$ for the incoherent state goes to zero according to $O(N−1/2)$ and we measured the $R$ for the system with size $N=1600$ for which $1/1600=0.025=O(10−2)$⁠, the chosen threshold value of $ε=0.02$ is appropriate.]. (ii) We estimated the critical value $pc=argmaxp(dR/dp)$ for which the slope of $R$ is maximized. We found that the values of $pc$ measured using both methods yield the same result and very well match our theoretical prediction Eq. (6) that we derived above, $pc=Q/(1+Q)$⁠, as is seen in the inset of Fig. 2.

The fully coherent state (⁠$R=1$⁠) is characterized by solutions on the coherence (or synchronization) manifold defined by $θi(t)=θj(t)=θs(t)$ for all oscillators $i$ and $j$⁠. Due to the phase-shift invariance, $θi(t)↦θi(t)+θ0$⁠, inherent to the model (1), there exists a reference frame co-rotating with the coherent solution in which $θi˙=0$ is satisfied for all oscillators $i$ and $j$⁠. Applying perturbations $εi(t)=θi(t)−θs(t)$ to the coherence manifold and linearizing Eq. (1), we obtain the variational equation

which is rewritten using the graph Laplacian $Lij=Kij−δijαi$ with Kronecker symbol $δij$ and row sum $αi:=∑j=1,j≠iNKij$⁠. The row sum of the Laplacian matrix is $∑j=1NLij=0$⁠, which implies the trivial eigenvector $(1,1,…,1)$ with eigenvalue $λ1=0$⁠; this is a consequence of the phase-shift invariance inherent to the system (1). Since the Laplacian matrix $L$ is real and symmetric, linear stability for the fully coherent (synchronized) state is guaranteed if all remaining $N−1$ (real) eigenvalues are less than zero (⁠$λi≤0$ for $i=2,3,…,N$⁠).

With this idea in mind, we computed the maximum eigenvalue $λmax<0$ of $L$ for a given value of $Q$ in order to determine the critical $pc$ for which linear stability of the coherent state breaks down; i.e., $λmax(pc)=0$⁠. Note that for the same value of $Q$⁠, many realizations for ${Kij}$ and thus also for ${Lij}$ can be chosen. We, therefore, computed an average value for $λmax$ obtained by ten distinct realizations of ${Lij}$⁠. Results for $λmax$ as a function of $p$ and varying values of $Q$ are shown in Fig. 3. We found that the agreement of these estimates for $pc$ with the theoretical prediction of $pc=Q/(1+Q)$ improves as the system size $N$ is increased.

In the following, we discuss an argument to determine the stability of the incoherent state along the line of Winfree’s and Kuramoto’s original work for all-to-all coupling with uniform strength (see Ref. 22 or 23, Secs. III and IV). Introducing the weighted order parameter (i.e., a mean-field),

where $Kij$ is given by (3), we can reformulate this field as follows: substituting (3) into (12) and using the definition for the Kuramoto order parameter (10) and $Q=−Kn/Kp$⁠, the field is now rewritten as

Note that we set $Kp=1$ and $Kn=−Q$ henceforth as mentioned further above. According to the mean-field theory on annealed networks,²⁴ we approximate the ER network as follows:

where $ki:=∑j=1NAERij$ is the degree of the $i$th node (or the $j$th node, respectively) and $⟨k⟩:=∑k=1NPkk$ is the mean degree with distribution $Pk$⁠. We may approximate the expectation value for the degree, $⟨k⟩≈N−1∑j=1Nkj=N−1Ne$⁠, where $Ne=N(N−1)p$ is the expected number of edges in $AER=AER(p)$⁠. For $N$ large, we can further approximate $ki≈kj≈⟨k⟩≈Np$ so that

(Note that this result coincides with the graphon for the ER network, as we discuss Sec. IV A further below.) Using this approximation in (13), we have

Finally, we formulate the governing equations (1) in terms of the mean weighted field $rieΦi$⁠, and we have $θ˙i=risin⁡(Φi−θi)$⁠; i.e.,

Expressing the phase dynamics in terms of the order parameter makes the mean-field character of the model obvious. We may now invoke a heuristic argument for the stability of incoherent oscillations.²³ Every oscillator appears to be decoupled from any other, although they all are interacting via the mean-field quantities $R$ and $Θ$⁠. Suppose there exists a quasi-stationary solution (incoherent or coherent) such that the order parameter is constant, $Rexp⁡(iΘ)=const.$ (this may be checked numerically or can be proven for the Kuramoto model with uniform coupling strengths). Then, an oscillator $i$ may be attracted to the angle of the mean-field rather than to the phase of any other oscillator; i.e., $|Θ(t)−θi(t)|→0$ as $t→∞$ if $(p(1+Q)−Q)>0$ but not otherwise.²⁵ Therefore, the critical value where the incoherent state loses stability is expected to occur at

Note that this value is identical to the critical value where the coherent state loses stability, Eq. (6). This suggests that hysteretic effects cannot be present, which is also confirmed by numerical computation.

Apart from disorder in the coupling, we also investigated how thermal noise in Eq. (1) affects the collective behavior of the system. Note that the two types of disorder influence the dynamics in different ways: while the disorder in the coupling, $Kij$⁠, remains constant in time (“quenched” disorder), the thermal noise [$ηi(t)$] models temporal fluctuations. For $T=0$⁠, we found that the system displays a discontinuous phase transition at a critical threshold $pc=Q/(1+Q)$⁠. How does thermal noise alter the characteristic of the phase transition, and what is the critical threshold for $T>0$?

To address these two questions, we performed numerical simulations of Eq. (1) and measured the order parameter $R$ in (10) while varying system parameters. We first varied the noise level $T>0$ while keeping the coupling ratio $Q=1$ fixed. This revealed features of a phase transition of second-order, as shown in Fig. 4. As temperature increases, so does the critical threshold $pc$ for the onset of synchrony increase. Next, we varied $p$ and $Q$ while keeping the temperature fixed at $T=0.1$⁠. The resulting behavior of $R$ is shown in Fig. 4. For the incoherent state (⁠$R=0$⁠) arising for $p<pc$⁠, we found that the averaged order parameter $R$ approaches zero as $N−1/2$⁠; for the partially coherent state arising for $p>pc$⁠, the order parameter $R$ approaches a saturated value $0<R<1$ as $N$ increases, such that any visible dependence on system size is not noticeable once the number of oscillators exceeds $N=200$⁠. The presence of thermal noise changes the nature of the transition from the incoherent (⁠$R=0$⁠) to the partially coherent state (⁠$R>0$⁠): for $T>0$⁠, the system exhibits a continuous transition between the two phases of second-order type, thus contrasting the first-order like transition found for the deterministic case with $T=0$⁠.

To numerically estimate the critical threshold $pc$ for this phase transition, we detected the value of $p$ for which the time-averaged order parameter exceeds the value of $R=0.1$⁠; i.e., we used the numerical threshold $ε=0.1$ to identify the incoherent state. (Note that the critical threshold $pc$ depends on the value $ε$⁠; i.e., a larger value of $pc$ results from larger values of $ε$⁠. Regardless, the critical threshold values $pc$ measured for $T>0$ exceed the values of those observed for $T=0$⁠).

The inset of Fig. 4 shows the behavior of $pc$ as a function of the coupling ratio $Q$ while keeping the temperature fixed at $T=0.1$⁠. We observed that the threshold $pc$ increased as temperature, $T$⁠, increased. This is plausible considering that larger thermal noise renders the synchronization of the system harder to achieve. Furthermore, increasing the value of coupling ratio $Q$ implies augmenting the dominance of the negative coupling. Thus, the system becomes harder to synchronize, and consequently, the critical threshold, $pc$⁠, becomes larger.

We apply recent results for a mean-field theory of the Kuramoto model valid for coupling networks with heterogeneous connectivity and coupling strengths based on the theory of graphons^1,26,27 or graphops.^2,28 The basic idea is that we may describe the phase evolution of oscillators in the mean-field limit, $N→∞$⁠, via a density $ρ(t,⋅)$ that evolves according to a transport equation, the Vlasov–Focker–Planck (VFPE) equation. This idea is originally formulated for complete graphs with uniform coupling strengths (full graphs); see, for instance, Refs. 23, 29, and 30.

We briefly outline how such a mean-field theory may be extended to describe heterogeneous, or even random, network topologies. To achieve this, it is instructive to consider how one can describe graphs in the mean-field limit, $N→∞$⁠. Consider sequences of weighted graphs $ΓN$ with increasing oscillator number $N$⁠, as illustrated in Fig. 5. These sequences may be thought of as an approximation of the description of the graph in the mean-field limit. Consider a discretization of $I:=[0,1]$ such that $Xn={ξ1(N),…,ξN(N)}$ and $limN→∞N−1∑i=1Nf(ξi(N))=∫If(x)dx$ for all continuous functions $f$ on $I$⁠. The points ${ξi(N)}i∈[N]$ are sampled either equidistantly, or randomly and uniformly i.i.d. Let $W$ be a symmetric and an integrable function on $I2=[0,1]×[0,1]$⁠. The weighted graph $ΓN=G(V,E)$ is defined by the node set $V(ΓN)=[N]$ and the edge set $E(ΓN)={(i,j)|W(ξi(N),ξj(N))≠0;i,j∈[N]}$⁠. Each edge $(i,j)$ is assigned the weight $Wij(N)=W(ξi(N),ξj(N))$⁠; if, on the other hand, the graph is random, we define a W-random graph via $P(i,j)=W(ξi(N),ξj(N))$ and 0 otherwise, where the decision for each pair $(i,j)$ is made independently. The limiting behavior of such sequences is determined by a symmetric measurable function on the unit square $I$ called a graphon. The graphon is a Fredholm integral operator defined by $W:L2(I)→L2(I)$⁠,

In the case of an ER graph (3), shown in Fig. 5, this sequence converges to a graphon with $WER(x,y)=p$⁠, as shown in Ref. 1. Note that this result represents the exact mean-field limit for $AER$ in Eq. (15).

Using graphons, it is possible to state a VFPE for a heterogeneous (random) graph structure in the mean-field limit. For simplicity, consider the deterministic case with a phase density $ρ(t,θ,x)$⁠. The density then evolves according to the transport equation,

with

with an initial condition $ρ(0,θ,x)=g(θ,x)$ sufficiently smooth. Here, $x$ can be thought of as a variable describing the heterogeneity of the network so that one effectively obtains a family of transport equations. One may take $x$ as a variable in the unit interval $x∈I=[0,1]$⁠, where points in the interval represent node labels in the infinite network limit.^1,27,31 It can be shown¹ that the density $ρ(t,⋅)$ is a unique weak solution for $t∈[0,T]$ of the initial value problem (20a)–(20b), characterized by an absolutely continuous measure $μt$⁠; similarly, that the solution of the initial value problem for the finite oscillators (1) generates a family of empirical measures $μtN$ for $t∈[0,T]$⁠; furthermore, these empirical and absolutely continuous measures converge as $N→∞$ for $t∈[0,T]$⁠, which implies convergence of the observed dynamics in the mean-field limit.

In Ref. 3, the authors develop a more general mean-field theory valid for the stochastic Kuramoto model; this means, in particular, that (20a) has an additional diffusive term that takes into account fluctuations of strength $T$⁠.³² The authors carried out a linear stability analysis for the incoherent solution of the stochastic Kuramoto model with $T>0$⁠. This is performed via linearizing the corresponding VFPE, resulting in an amplitude equation describing the growth (or decay) of perturbations around the incoherent branch. For the special case of the standard Kuramoto model with uniform coupling, $Kij=K$ for all $i,j∈[N]$⁠, and identical frequencies, one obtains the critical coupling for the incoherence–coherence transition,

a result that coincides with the classical result by Sakaguchi³³ and later by Strogatz and Mirollo.³⁴ 

The coupling in (3) can be expressed in terms of ER networks—accordingly, the corresponding graphon valid for the mean-field limit (⁠$N→∞$⁠) is also a linear combination of a constant and the graphon $WER$ for the ER network,

The last equation follows as the graphon for the ER network has been derived in Ref. 1 to be $WER(x,y)=p$⁠. Note that in this case, random sampling of the (finite) graph $ΓN$ must be done by modifying the prescription for the W-random graph introduced further above; i.e., each edge $(i,j)$ is attributed to the weight via $P(i,j)=Wp,Q(ξi(N),ξj(N))$ and $W0,Q$ otherwise.

Since we without loss of generality chose $Kp=1$ (and $Kn=−Q$⁠) and the graphon (22) is a constant that can be placed outside the integral in (20b), we effectively have an effective coupling strength $Keff=−Q+(1+Q)p$⁠. Letting this effective coupling be critical, $Keff=Kc$⁠, we obtain from (21) a condition for the change of stability of the incoherent branch,

The simulations and measurements of the threshold for the incoherence–coherence transition shown in Fig. 4 show an excellent agreement with (23). Note that our previous result in (18) is retrieved as $T→0$⁠. Finally, the same result for the critical threshold could also be obtained by calculation of the eigenvalues of the graphon;^1,3 however, the above solution is more direct.

We considered a population of all-to-all coupled oscillators subject to quenched and temporal disorder and explored its collective behavior. The quenched disorder arises due to (symmetric) random coupling interactions with either a positive or negative value and the temporal disorder due to additive thermal noise. For the deterministic case (⁠$T=0$⁠) with quenched disorder only, we found that the system displays a discontinuous first-order like phase transition from the incoherent to the fully coherent state. We provided theoretical arguments for the nature of this transition and explained its critical threshold value $pc$ for the thermodynamic limit. We confirmed the critical threshold $pc$ rigorously by providing a stability argument and performing a linear stability analysis of the fully coherent state (⁠$R=1$⁠). Vice versa, for the noisy case (⁠$T>0$⁠) with additional temporal disorder, we found that the system displays a continuous second-order phase transition from the incoherent to the partially coherent state (⁠$0<R<1$⁠). Numerically, we found that the corresponding critical threshold $pc$ grows with increasing $Q$ and $T$⁠, and we determined an exact formula for the incoherence–coherence transition based on a rigorous mean-field theory³ and a linear stability analysis for the incoherent branch. The formula for this threshold reduces to the results we obtained for $T=0$⁠.

A number of studies concerned the collective dynamics in networks of coupled oscillators characterized by heterogeneous properties and interactions. Such models may have contiguous properties such that oscillators form subpopulations or communities with identical properties, giving rise to multimodal frequency distributions or non-uniform coupling interactions demarcating subpopulations of strongly interacting oscillators.^29,35–39 Another type of heterogeneity arises when these properties are non-contiguous, including a variety of oscillator models with heterogeneous positive and negative coupling. For instance, our model is related to the important class of oscillator glass models,^8–13 which distinguish themselves from (1) in that they assume distributed natural frequencies and normally distributed coupling strengths $Kij$⁠. A special case represents the $XY$ spin-glass models for identical oscillators with $ωi=0$⁠; see Ref. 7. To render these problems analytically tractable, most studies assumed separable disorder in the coupling; e.g., $Kij=Jsisj$ with $si(j)=±1$ and $J$ a constant. In this case, the spin-glass phase is found to exist in the regime of the incoherent state where the $XY$ spins are fully randomized.⁷ For a review on such types of models, see also Ref. 40. A simpler coupling scheme arises in systems with oscillator nodes that have uniform input weights; i.e., $Kij=Ki$ for all $i∈[N]$⁠.^41,42 Hong and Strogatz^15,16 studied this case where coupling strengths are either value $Ki=Kp>$ or $Ki=Kn<0$⁠. The oscillators were found to undergo a phase transition from the incoherent (⁠$R=0$⁠) to the partially coherent state (⁠$0<R<1$⁠) when natural frequencies were normally distributed; vice versa, if natural frequencies were identical, oscillators were found to display a phase transition from an incoherent state to a state exhibiting a traveling wave.¹⁵ Alternatively, one may consider uniform output weights (⁠$Kij=Kj$ for all $j∈[N]$⁠); for this case, it was found that in the presence of heterogeneous frequencies, the system displays conventional mean-field behavior with a second-order phase transition like that found in the original Kuramoto model.¹⁷ An intriguing extension of coupled oscillator models is the swarmalator model,⁴³ where oscillators are allowed to move on a spatial domain, where their respective level of synchronization influences their motion, while their distance modulates their respective coupling strengths. The collective behavior resulting from such a model with the same type of quenched random disorder studied here has recently been investigated.⁴⁴ 

Remarkably, the critical threshold $pc$ [Eq. (6)] for quenched disorder (⁠$T=0$⁠) is identical to the one observed for the systems mentioned above, i.e., including uniform input weights¹⁵ (see also Ref. 45 with $Q=1$ and system parameter $γ=0$ therein), uniform output weights,¹⁷ as well as the model of identical swarmalators⁴⁴ with quenched random disorder. Indeed, upon further inspection, the universality of the threshold $pc$ for $T=0$ seems fairly evident from the energetic argument we made. Specifically, the energy for the completely incoherent state is $H=0$ in the thermodynamic limit $N→∞$⁠, while it is negative for the fully coherent state, $H<0$⁠. Assuming that the energy density changes continuously, we may identify the critical state as the point where the energy density of the two phases becomes equal. As we have seen, this critical transition is closely related with the mean value of the coupling disorder, given by

the value of which also crosses zero at criticality since $Ecoh=−⟨Kij⟩$⁠. This argument yields the critical threshold $pc=Q/(1+Q)$ and does not depend on the specific structure of the coupling matrix $Kij$ (i.e., the ordering in terms of the values $Kn$ and $Kp$⁠), but rather only on the specific number of entries with negative or positive coupling. Note that the general coupling matrix $Kij$ includes all cases discussed here, including uniform input (⁠$Kij=Ki$⁠) or output (⁠$Kij=Kj$⁠). Thus, this explains the universal nature of the critical threshold for the various quenched types of random disorder in the coupling strength discussed here. However, while this universality of the critical threshold holds when thermal noise is absent (⁠$T=0$⁠), our numerical simulations indicated that this universality breaks down when thermal noise is present (⁠$T>0$⁠).

While most classical studies assumed separable random disorder for the coupling interactions, we considered non-separable coupling interactions. For future research, it would be interesting to investigate whether the critical threshold $pc$ for the noisy case (⁠$T>0$⁠) also has a universal character as the zero-noise case (⁠$T=0$⁠) in terms of various types of temporal disorder in the coupling. We leave these question for a future study.

We would like to thank Hyunggyu Park and Kangmo Yeo for useful discussions in the early stage of the study, Jaegon Um and Christian Kuehn for helpful discussions regarding the stability analysis, and Christian Bick regarding valuable comments on the manuscript. This research was supported by the NRF (Grant No. 2021R1A2B5B01001951) (H.H.).

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.