Validity of annealed approximation in a high-dimensional system

This study investigates the suitability of the annealed approximation in high-dimensional systems characterized by dense networks with quenched link disorder, employing models of coupled oscillators. We demonstrate that dynamic equations governing dense-network systems converge to those of the complete-graph version in the thermodynamic limit, where link disorder fluctuations vanish entirely. Consequently, the annealed-network systems, where fluctuations are attenuated, also exhibit the same dynamic behavior in the thermodynamic limit. However, a significant discrepancy arises in the incoherent (disordered) phase wherein the finite-size behavior becomes critical in determining the steady-state pattern. To explicitly elucidate this discrepancy, we focus on identical oscillators subject to competitive attractive and repulsive couplings. In the incoherent phase of dense networks, we observe the manifestation of random irregular states. In contrast, the annealed approximation yields a symmetric (regular) incoherent state where two oppositely coherent clusters of oscillators coexist, accompanied by the vanishing order parameter. Our findings imply that the annealed approximation should be employed with caution even in dense-network systems, particularly in the disordered phase.


Introduction
Recently, there has been notable attention given to dynamics of complex systems.One popular strategy for understanding their connection geometry is through the use of networks composed of nodes and links 1,2 .Network structure often exhibits quenched link disorder, thereby rendering the system analytically intractable.In addition to numerical analyses, mean-field approximations have been employed to study the collective properties of complex-network systems.The "annealed" approximation (AA), a frequently employed mean-field method [1][2][3][4][5][6][7][8][9][10][11][12][13][14] , characterizes a quenched link as an annealed one with an appropriate linking probability.This approximation is commonly referred to the heterogeneous mean-field approximation [1][2][3] in complex-network studies, where the linking probability depends solely on the numbers of links (degrees) of the two connecting nodes.
In sparse networks characterized by a finite mean degree, it is well-established that the AA manifests various limitations [11][12][13][14][15] , as it only captures a portion of network disorder.For instance, in investigations of Kuramoto-type models of coupled oscillators [16][17][18][19] , previous studies have observed that systems subject to the AA not only shift the transition point but also occasionally alter the nature of the transition 14 .On the contrary, in the context of dense networks with a diverging mean degree 20 , there exists a widely held belief that such networks closely approximate a complete graph (all-to-all connections) in the thermodynamic limit.Consequently, it is reasonable to anticipate that the AA may correctly characterize the collective properties of dense-network systems in general, as the AA mitigates disorder fluctuations in such systems 21 .It is reminiscent of the adiabatic elimination method commonly employed in quantum optics 22 : According to the adiabatic elimination, the dense-network fluctuations can be regarded as fast modes, resulting in corrections to the slow-mode dynamics corresponding to the fully connected case.Nevertheless, in the incoherent (disordered) phase where the order parameter vanishes, finite-size effects may exert significant influence on stabilizing steady states, suggesting a breakdown of the adiabatic elimination because of a lack of slow-modes.Thus, one might imagine a potential disparity in the incoherent steady-state patterns for systems on the complete graph (CG), dense network (DN), and annealed network (AN), attributed to the finite-size effects contingent upon the network structure.
In this study, we explore the extent of similarity between a DN system and its CG counterpart, and assess the validity of the AA.To achieve this, we employ systems of the Kuramoto-type oscillators known for exhibiting collective properties sensitive to connectivity (link) disorder, which has led to a critical failure of the AA in sparse-network systems 14 .Firstly, we demonstrate that connectivity fluctuations in DNs vanish in the thermodynamic limit, resulting in the order parameter behavior identical to that of the CG version.In the case of identical oscillators without frequency disorder, often referred as the arXiv:2403.14984v1[cond-mat.stat-mech]22 Mar 2024 Watanabe-Strogatz (WS) model 23,24 , we observe that finite-size effects stemming from connectivity disorder are strong enough to destroy all infinitely many, initial-condition dependent, "regular" (symmetric) steady-state patterns found in the incoherent phase of its CG version 23,24 .In DNs, random steady-state patterns emerge in the end, suggesting that quenched disorder can readily eliminate regularity in incoherent patterns.Under the influence of the AA, similar incoherent regular patterns appear with more complex symmetries, distinct from both the random patterns observed in DNs and the regular patterns in the CG version.We note that all these regular patterns are vulnerable to frequency disorder, resulting in random incoherent patterns.Consequently, the incoherent patterns are all identical in the CG, DN, and AN with heterogeneous oscillators featuring random natural frequencies.
Secondly, we introduce coupling disorder (competition of attractive and repulsive couplings) to the WS model on the CG and examine its collective behavior both numerically and by applying the AA.This model can be also viewed as the WS model on a combined version of two DNs based on two competing couplings.Intriguingly, we discover a single initial-condition independent symmetric steady-state pattern in the incoherent phase for the system subject to the AA.In this symmetric pattern, oscillators are sharply divided into two coherent clusters with a phase difference of π.The sizes of these two clusters become identical in the thermodynamic limit, resulting in the vanishing order parameter.Our numerical analysis without the AA reveals that this incoherent symmetric pattern is replaced by the incoherent random pattern.This sharp and simple disparity in incoherent states leads us to conclude that the AA may not accurately represent the incoherent steady-state pattern even in (generalized) DNs for a wide range of various many-body dynamics.
Finally, we note that temporal networks where connection links undergo temporal changes bear resemblance to the annealed network systems when the time scale for link connection is sufficiently short.Thus, in these networks which are prevalent in various biological and social systems, the regular/symmetric incoherent patterns can be empirically observed.

Kuramoto model on dense networks
We start with a system of N Kuramoto-type oscillators 16 on a complex network.The dynamics of each oscillator is governed by the equation of motion as where φ i and ω i represent the phase angle and natural random frequency of oscillator i, respectively.The element a i j of the adjacency matrix denotes the connectivity between oscillators i and j, with a i j = 1 indicating a connection and a i j = 0 otherwise.The parameter J represents the strength of coupling.For simplicity, the summation is normalized by the mean degree ⟨k⟩ = (1/N) ∑ i k i , where k i (= ∑ j a i j ) denotes the degree (the number of neighbors) of oscillator i.
Defining the "local" field h i by with phase factor z j = e iφ j , Eq. ( 1) is now rewritten as where Im(X) denotes the imaginary part of X and z * i the complex conjugate of z i .We also define the order parameters from averaged z i and h i , respectively: Note that ⟨z⟩ is the so-called Kuramoto phase order parameter defined in Ref. 16.
It is convenient to recast the local field in the form as where δ ki and ξ i represent the degree fluctuation and the average phase factor (z i ) fluctuation over neighbors, respectively, with and a "noise" term ξ i reads In the case of the CG where a i j = 1 for all pairs, the local field and two order parameters become identical, i.e., h i = ⟨h⟩ = ⟨z⟩ with the complete absence of fluctuations (δ ki = 0 and ξ i = 0).In contrast, for a sparse network characterized by a finite ⟨k⟩ such as the Erdös-Rényi (ER) network with an extremely low connection probability 25,26 or scale-free networks 27 , the fluctuations can be of the order O(1).Consequently, these fluctuations may influence the dynamics significantly, leading to steady states distinct from those of the CG version.
Here, we focus our attention on a dense network (DN) characterized by the following conditions: with 0 < α ≤ 1 and β < α.The above power-law scaling of ⟨k⟩ can be found in the ER network with a finite connection probability 26 and scale-free-like DNs 27 .It is then straightforward to estimate the order of local fluctuations as where we estimate ξ i ∼ O( √ k i /⟨k⟩) under the reasonable assumption that the phase factor fluctuation (z j − ⟨z⟩) is almost Gaussian random with zero mean.In DNs, both local fluctuations vanish as N → ∞, and then the local field h i in Eq. ( 5) converges to the global order parameter ⟨z⟩ in the ordered phase with ⟨z⟩ ∼ O(1).Consequently, the dynamic equation of motion, Eq. ( 3), becomes identical to that for the CG as N → ∞ and thus the order parameter values for DNs become identical to those for the CG.It is noteworthy that in the incoherent (disordered) phase, where ⟨z⟩ also approaches to zero as N → ∞, ⟨z⟩ competes with the local fluctuation ξ i in Eq. ( 5).This suggests that different types of finite-size effects in the CG, DN, and AN may yield distinct incoherent patterns, while still adhering to the condition of the vanishing order parameter.

Annealed approximation on networks
The analytical treatment of models on a network with quenched link disorder is typically challenging.Instead, the AA is frequently employed due to its analytical tractability, while still yielding results analogous to the original system.
In the AA, networks are substituted with annealed ones using the heterogeneous mean field approximation [1][2][3] , where the linking probability depends solely on the degrees of the two connecting nodes.Within this framework, the adjacency matrix a i j is replaced by a A i j as which represents the mean linking probability between nodes i and j 3,6 .Then, the annealed version of the local field h A i is given by with phase factor z A j = e iφ A j satisfying the AA-applied dynamic equation.The average phase factor fluctuation ξ A i in the AA is rather simplified as Note that the local property of ξ A i is limited to k i (no explicit connection information), leading to the simple expression for h A i as whose local property is also solely given by k i , allowing the analysis of the dynamic equation ( 3) simpler.The order of ξ A i is estimated similarly as which also vanishes as N → ∞ for DNs, as expected.Nonetheless, the finite-size effects are different for the annealed and quenched networks, which may exert a crucial influence on the determination of incoherent steady states.

Validity of the annealed approximation for the WS model
For the WS model with identical oscillators (ω i = Ω) 23,24 , the dynamic equation is given by φi = Ω + J Im (h i z * i ).We can set Ω = 0 without loss of generality via simple mapping of φ i → φ i + Ωt.On the CG, the local field loses its locality completely as In the long time limit, there exist two types of stable fixed points {φ i }, satisfying (a) ⟨z⟩ = 0 for J < 0 (see Methods) or (b) φ i = Φ for J > 0 with the global angle Φ defined by ⟨z⟩ = |⟨z⟩|e iΦ .Consequently, we obtain |⟨z⟩| = 1 for J > 0 and ⟨z⟩ = 0 for J < 0 with {φ i } satisfying discrete rotational symmetries (∑ N j=1 z j = 0) 23,24 .An initial condition selects one of these many incoherent "regular" symmetric steady states for J < 0, exact for any finite N.
In DNs, these incoherent regular fixed points are not stable any longer, due to strong fluctuations in the local field h i .Instead, φ j becomes random with the uniform distribution over [0, 2π] to make ⟨z⟩ vanish in the N → ∞ limit for J < 0. This is confirmed by numerical simulations [not shown here].In the AA of DNs, the local field fluctuations become weaker with a degree dependence only as in Eq. ( 13), leading to Similar to the CG case, there are two types of stable fixed points, satisfying (a) ⟨h A ⟩ = 0 for J < 0 (see Methods) and (b) φ A i = Ψ A for J > 0 with another global angle Ψ A defined by ⟨h A ⟩ = |⟨h A ⟩|e iΨ A .Consequently, we obtain |⟨z A ⟩| = 1 for J > 0 and ⟨z A ⟩ → 0 as N → ∞ for J < 0 with {φ A i } satisfying more complex discrete symmetries (∑ N j=1 k j z A j = 0).Again, the incoherent steady state depends on the initial condition.
The incoherent steady state patterns are all different for the WS model on the CG, DN, and AN.It implies that the AA fails to predict the stable steady-state solutions of Eq. (1) for DN systems, in particular when the oscillators are identical.Certainly, this failure is not common in general DN systems.An illustrative example is evident in the Kuramoto model with frequency disorder in ω i .In its incoherent phase in DNs, oscillators undergo rotational motion with their respective natural frequencies, i.e. φi = ω i in the long time limit.The order parameter ⟨z⟩ as well as the local field h i vanish as N → ∞, identical to the behavior observed in the CG.The AA-applied systems are expected to yield the same rotating behavior, thereby resulting in identical incoherent states on the CG, DN, and AN.This is not much surprising, as the additional frequency disorder is of the order O(1), dominant over local fluctuations that diminish in the N → ∞ limit for DN systems.Temporal fluctuations such as a thermal noise also destroy the intricate structure of incoherent regular steady states in the AA (and also in the CG), and then the incoherent steady state pattern becomes fully random, akin to the incoherent phase in DN systems.
An intriguing question arises as to the nature of additional disorder that maintains a distinction in incoherent steady-state patterns between the AA-applied and original system.Subsequently, we investigate the inclusion of a quenched coupling disorder residing in the links of DNs, where the associated fluctuations are expected to be of a similar order to those of the quenched link disorder.

Oscillators with competing interactions
We consider a generalized Kuramoto model with coupling disorder, governed by the dynamic equation as where J i j is a random coupling strength between oscillators i and j.For simplicity, we consider the CG case (a i j = 1, thus ⟨k⟩ = N).This model has been introduced and studied in the context of "oscillator glass" by Daido 28,29 .
In the case of identical oscillators (ω i = 0), this model describes the zero-temperature Sherrington-Kirkpatrick model for XY spin glass [30][31][32] .Recently, Hong and Martens 33 also considered this model with a probability distribution P[J] for J i j that has two δ -peaks such as with J + > 0 (attractive) and J − < 0 (repulsive).These competitive interactions induce frustration between oscillators and the first-order phase transition occurs between the fully ordered and disordered phases when the mean value of coupling strengths (pJ + + (1 − p)J − ) changes its sign, i.e. at p = p c = (−J − )/(J + − J − ) 33 .
It is now worth noting that the random interaction J i j in the so-called "two-peak" model with identical oscillators can be expressed using the adjacency matrix b i j of a random network, as follows: where b i j = 1 represents a positive (attractive) link with J + , while b i j = 0 represents a negative (repulsive) link with J − .As the positive and negative links are randomly distributed, the positive degree k i (the number of positive links stemming from node i) satisfies the binomial distribution with mean ⟨k⟩ = N p and variance ⟨k 2 ⟩ − ⟨k⟩ 2 = N p(1 − p) for large systems.Thus, the {b i j } network is dense with the exponents α = 1 and β = 1/2 defined in Eq. ( 8).This model can be also regarded as a competition between the CG model with the negative coupling constant J − (< 0) and the DN model (b i j ) with the positive one p(J + − J − ) (> 0).For more general case with the underlying DN (a i j ) and multi-peak distributions of P[J] in Eq. ( 17), we expect the coupled systems on multiple DNs like a hypernetwork 34 , which are left for future study.Eq. ( 19) is rewritten in a more illustrative form as φi where the corresponding local field q i reads with ∆ ≡ (−J − )/ (J + − J − ) = p c (0 < ∆ < 1).Using Eq. (5), we get with ξ i = ∑ j b i j (z j − ⟨z⟩) /⟨k⟩.As discussed previously, both local fluctuations, δ ki and ξ i , in the DN diminish as N → ∞, thus the local field q i approaches ⟨z⟩(p − ∆).Consequently, the dynamic equation becomes identical to Eq. ( 15) with substitution of J by (J + − J − )(p − ∆) = pJ + + (1 − p)J − (mean coupling strength).Therefore, we expect the first-order phase transition at p = ∆ from the disordered to the fully ordered phase.As the {b i j } network is dense (not CG for p < 1), there are finite-size fluctuations which make the incoherent steady state uniformly random for p ≤ ∆ (numerically confirmed later, shown in Fig. 1).
The application of the AA as in Eq. ( 10) yields with the annealed local field q A i given by where Eq. ( 13) is used.Similar to the previous cases without coupling disorder shown in Eqs. ( 15) and ( 16), there exist incoherent regular fixed points satisfying both ⟨h A ⟩ = ⟨z A ⟩ = 0 simultaneously.However, these regular fixed points are proven to be always unstable with a distribution of k i (see Methods).Instead, one of two "ordered" fixed points becomes stable in the long-time limit, depending on the magnitude of degree k i ; with ⟨z A ⟩ = |⟨z A ⟩|e iΦ A .The stability condition for the fixed points requires the identical global phase angles, denoted as Φ A = Ψ A where ⟨h A ⟩ = |⟨h A ⟩|e iΨ A .Thus, the threshold value k * is determined by the equation k * |⟨h A ⟩|/N − ∆|⟨z A ⟩| = 0 (see Eq. ( 24)), which can be solved in a self-consistent manner.Detailed derivations are given later.

5/11
The coexistence of two coherent (ordered) clusters with a phase difference of π is observed across the entire parameter range of (p, ∆).In the disordered phase (p < ∆), the threshold value k * ≃ ⟨k⟩, leading to identical cluster sizes in the N → ∞ limit.Thus, the order parameters approach zero; ⟨z A ⟩ ≃ 0 and ⟨h A ⟩ ≃ 0. Conversely, in the ordered phase (p > ∆), one cluster fully dominates over the other in the N → ∞ limit with k * ≃ (∆/p)⟨k⟩.This results in ⟨z A ⟩ ≃ 1 and ⟨h A ⟩ ≃ 1.
It is evident that the incoherent steady state for the AA-applied system differs from the random incoherent one in the original system.This incoherent steady-state pattern exhibits the simple Z 2 symmetry (the same number of oscillators with z A i e −iΦ A = 1 or −1), resembling a regular symmetric pattern found in the WS model on the CG in Eq.( 15).However, the origin of this symmetry is obviously different.Furthermore, this symmetry is exact only in the N → ∞ limit and is also independent of initial conditions, while the regular symmetric patterns in the previous models are exact for any finite N and are dependent on initial conditions.i and yellow ones represent φ i .The phase angle segregation in terms of degrees emerges in the annealed system, while it does not in the quenched case.Note that the mean degree is given by ⟨k⟩ = N p = 640.

Numerical results
We numerically solve the equation of motion, Eq. ( 19), employing the Heun's method with various network sizes and fixed p = ⟨k⟩/N=0.2.Initial phase angles are randomly chosen within [−0.005π, 0.005π].Solid symbols in Fig. 1 (a) denote the values of |⟨z⟩|, averaged over a period of t = 5 × 10 4 ∼ 10 5 , after discarding an initial transient period of the same duration and also averaged over 10 ∼ 100 network realizations and initial conditions.We also display the analytic solutions |⟨z A ⟩| in Fig. 1 (a), obtained from a self-consistency equation Eq. ( 35) for the annealed networks.We find that both |⟨z⟩| and |⟨z A ⟩| seem to show a discontinuous transition from the fully ordered to the disordered phase at ∆ = p, as |⟨z⟩| tends to approach 1 (0) asymptotically for ∆ < p (∆ > p) for increasing N.At ∆ = p, however, |⟨z A ⟩| remains finite as N → ∞, as demonstrated by the crossing of solid curves, while |⟨z⟩| approaches zero.This observation may suggest distinct underlying mechanisms of the transitions between the annealed and quenched systems.
In (b)-(e) in Fig. 1, phase angle snapshots are plotted as a function of degree for the quenched and annealed cases, denoted by φ and φ A , respectively.Data are obtained numerically from Eq. ( 19) in a single network realization for p = 0.2 and N = 3200.Starting from an initial condition described above, the data points are obtained at t = 10 5 near the steady state.For ∆ < p, we observe in Fig. 1 (b) and (c) that in the original quenched system a single coherent cluster with φ ≃ 0 is formed by oscillators with higher degrees (stronger interactions effectively) and scattered phase angles for those with lower degrees seem to be due to finite-size effects.The annealed case shows two coherent clusters with phase difference of π, but one cluster dominates over the other.The finite-size effects are much weaker in the annealed systems, as their fluctuations should be much weaker than those for the quenched systems.
For ∆ > p in Fig. 1 (e), the phase angles seem randomly distributed as expected in the quenched system.Remarkable distinction is found in the annealed system, where a binary mixture of two coherent clusters with comparable sizes emerges.The contributions of the two clusters to the order parameter ⟨z A ⟩ seem to cancel out exactly in the N → ∞ limit.At the transition (∆ = p), the balance of two cluster sizes are slightly broken, yielding a nontrivial value of |⟨z A ⟩| ≃ 0.61946, which is consistent with the analytic result derived in the following.
Back to the scattered plots in Fig. 1, we note that the outcomes from the original quenched model (orange dots) present notable disparities in comparison to those of the annealed model (blue dots), particularly in the vicinity of the threshold value k * .This phenomenon can be understood by considering the following.In the quenched system, the magnitude of noise fluctuations ξ i is of the order O(N −1/2 ).This suggests that the local field q i is predominantly influenced by these fluctuations, especially near k ≃ k * , where the value of the annealed local field almost vanishes as q A i ≃ 0. In contrast, in the region where ), the annealed local field q A i becomes comparable to ξ i , resulting in a reduced disparity between the quenched and annealed results.This trend is also observable in the scattered plots.

Analytic solutions of the annealed two-peak model
We rewrite the dynamic equation ( 23) in a convenient form as where yielding the stable steady-state fixed points as Utilizing the definition of the global phase angle Φ A , we obtain the expression as |⟨z A ⟩| = z A e −iΦ A = ⟨e i φ A ⟩, implying ⟨sin φ A ⟩ = 0 and ⟨cos φ A ⟩ ≥ 0. In addition, from Eqs.( 27) and ( 28), we observe that the signs of the steady-state solution for sin φ A i should coincide with the sign of sin(Ψ A − Φ A ), independent of i.Then, the constraint of ⟨sin φ A ⟩ = 0 demands sin φ A i = 0 for all i's, thereby A i = 0 and sin(Ψ A − Φ A ) = 0. Subsequently, it follows that cos φ A i = ±1 and cos(Ψ A − Φ A ) = ±1 in the steady state.With the selection of cos(Ψ A − Φ A ) = −1, B i is always negative, thus Eq. (28) yields cos φ A i = −1 for any i.However, this contradicts the condition ⟨cos φ A ⟩ ≥ 0, thereby cos(Ψ A − Φ A ) = +1 (implying Ψ A = Φ A ) must be selected for stable fixed points.In this case, the sign of B i changes as k i varies.Consequently, cos φ A i = +1 with B i > 0 for large k i and cos φ A i = −1 with B i < 0 for small k i , leading to where k * is defined as To determine the threshold value k * , we need to solve Eq. ( 30) in a self-consistent manner.Using |⟨z A ⟩| = ⟨cos φ A ⟩ and |⟨h A ⟩| = ⟨k cos φ A ⟩/⟨k⟩, the order parameters in the steady state are given as

7/11
where the degree distribution B(k, p) is the binomial distribution of degree k for a given p.For large N, it is well known that B(k, p) can be approximated by the Gaussian distribution of a continuous variable k ∈ [−∞, ∞] with mean ⟨k⟩ = N p and variance N p(1 − p).In this continuum limit, Eq. ( 31) can be expressed in a simple form as where the error function is defined as erf[x] = 2/ √ π x 0 dt exp −t 2 .From the constraint |⟨z A ⟩| ≥ 0, we note that k * cannot exceed ⟨k⟩ = N p.
For convenience, we rewrite Eq. ( 32) as with By substituting Eq.( 33) into Eq.( 30), we obtain erf The solution for X in Eq. ( 35) provides the threshold value k * and also determine the order parameter values in Eq. ( 33) in the steady state.These results are in very good agreement with the numerical outcomes obtained directly from the annealed equation ( 23) (not shown here).For small ε (large N), the explicit expression for X can be derived from Eq. ( 35) as with X c ≈ −0.62006, determined by solving the equation of erf [X c ] e X 2 c X c = 1/ √ π which is given by Eq. ( 35) with ∆ = p for small ε.
In large quenched systems, the numerical results suggest ρ (φ ) = δ (φ − Φ) for ∆ < p.In contrast, for ∆ ≥ p, it is observed that the phase angles are uniformly distributed, implying that ρ(φ ) = 1/(2π).We can also consider the generalized order parameters defined as z m ≡ |⟨z m ⟩| and z A m ≡ |⟨ z A m ⟩| with an arbitrary integer m.We note that z A m for even m does not distinguish the ordered and disordered phase (z A m = 1) in the annealed systems, as two coherent clusters are synchronized with phase difference of π.

Summary
In summary, we showed that the dynamic equations governing DN systems with quenched link disorder converge to those of the CG version in the thermodynamic limit, where the local fluctuations vanish entirely.Consequently, the AA-applied systems where fluctuations are attenuated, exhibit the same dynamic behavior in the thermodynamic limit.However, a notable discrepancy may arise in the incoherent (disordered) phase, where finite-size effects can become critical in determining the steady-state pattern.We illustrate our findings through two prototypical models of coupled oscillators; the WS model for synchronization and the zero-temperature XY-like model with competing couplings.In both cases, we analytically derive the incoherent patterns in the annealed systems, revealing stark differences from those in the original quenched systems.These patterns in the annealed systems are regular and symmetric, in contrast to the random patterns in the quenched case.This suggests that caution should be given when applying the AA even to the DN systems, particularly when examining incoherent steady-state patterns.We emphasize that our analysis based on the network topology and findings are not restricted to the oscillator or XY-spin systems.One may expect a similar discrepancy for the two-peak model with Ising spins because phase angles of oscillators with the AA-applied two competing interactions align along two branches, reminiscent of the Ising spins.Furthermore, we point out the possibility of observing these intriguing regular and symmetric incoherent patterns in temporal networks with sufficiently fast time scales, which may underlie various biological and social systems.

Stability analysis of incoherent regular solutions
First, we perform a linear stability analysis of the incoherent regular fixed points for the WS model.For the CG version with the dynamics governed by Eq. ( 15), a small perturbation δ φ i around the fixed points φ s i satisfying ⟨z⟩ = 0 evolves as with the stability matrix element S i j given by This stability matrix can be expressed as the sum of two rank-1 matrices as S i j = S (1) i j with S i j = (J/N) cos(φ s j ) cos(φ s i ) and S (2) i j = (J/N) sin(φ s j ) sin(φ s i ).It is trivial to show that any rank-1 matrix has zero eigenvalues except for one eigenvalue given by its trace.For J > 0, both rank-1 matrices are positive semi-definite with the non-zero eigenvalues of (J/N) ∑ j cos 2 (φ s j ) and (J/N) ∑ j sin 2 (φ s j ), respectively.As the sum of two positive semi-definite matrices should be also positive semi-definite, it follows that the eigenvalues of the stability matrix are non-negative, leading to the conclusion that all regular fixed points should be unstable for J > 0. Consequently, the stabilization of the regular fixed points is only feasible for J < 0.
This stability analysis can be also applied to the AN version with the dynamics governed by Eq. ( 16), leading to the stability matrix S A i j as It is clear that this stability matrix can be also decomposed into two rank-1 matrices, thus the same conclusion as above can be drawn, as the link degree k i is always non-negative.
For the AA-applied model with coupling disorder, the dynamic equation is given by Eq. ( 23), yielding the stability matrix as It should be noted that this matrix can be decomposed into two positive semi-definite rank-1 matrices and two negative semi-definite rank-1 matrices.It is straightforward but rather lengthy to prove that the signs of the non-zero eigenvalues of the stability matrix align with those of the decomposed rank-1 matrices, if none of these rank-1 matrices is not a multiple of another (a proof not included here).Given a distribution of k i , it is evident that none of our four rank-1 matrices is proportional to any other.Consequently, the stability matrix invariably possesses two positive eigenvalues, thereby rendering all incoherent regular solutions unstable.This conclusion is further substantiated by numerical results.

Figure 1 .
Figure1.The original quenched system of the two-peak model versus the annealed-network version.Data points were obtained from numerical integrations of Eq. (19) with the quenched network b i j and its annealed version for p = 0.2 and J + = 1.(a) The order parameter as a function of ∆ plotted on a semi-log scale for various values of N. Solid symbols denote |⟨z⟩| for the original system and the dashed lines are guides to the eyes.Solid lines represent the analytic solutions of |⟨z A ⟩| obtained by the self-consistency equation (35) for annealed networks.Both the analytic and numerical results suggest a discontinuous transition at ∆ = p as N → ∞.The order parameter |⟨z⟩| at the transition seems to approach zero, while its annealed version |⟨z A ⟩| converges to a nontrivial value.(b)-(e) Snapshots of phase angles φ and φ A versus degree k in a single network realization of N = 3200 near the steady state for various values of ∆.Blue dots represent φ A i and yellow ones represent φ i .The phase angle segregation in terms of degrees emerges in the annealed system, while it does not in the quenched case.Note that the mean degree is given by ⟨k⟩ = N p = 640.

Figure 2 .
Figure 2. Schematic diagram summarizing phase-angle distributions in annealed and quenched systems.Black solid lines represent the generalized order parameters z m = |⟨z m ⟩| for the quenched system with integer m.Blue solid lines represent the generalized order parameter z A m for the annealed systems with odd m, while red solid lines with even m.