A Parameterized locked generator phases

# Mean-field approach for frequency synchronization in complex networks of two oscillator types

## Abstract

Oscillator networks with an asymmetric bipolar distribution of natural frequencies are used to describe stable operations in electrical transmission grids. We propose a mean-field model that captures the onset and form of frequency synchronization in such oscillator networks. The model takes into account a broad class of heterogeneous connection structures. It yields critical thresholds and a closed-form description of the synchronized state, identifying basic properties that synchronized regimes possess on all considered topologies. The framework also captures synchronized regimes with large phase differences that commonly appear just above the critical threshold. Additionally, two model quantities gauge the accuracy of mean-field assumptions.

## I Introduction

Synchronization is a frequent phenomenon in nature, appearing in biological, ecological, sociological, and engineering contexts (see [1]; [2]; [3]; [4]; [5] for surveys). These systems can be considered networks of coupled phase oscillators and are described by the paradigmatic Kuramoto model [6] or its extensions. In this setting, the phase of oscillator is driven both by its natural frequency and by the phase difference with adjacent oscillators that couple with uniform strength . The connection structure is encoded in the adjacency matrix with entries if oscillators and are connected and otherwise, with the dynamics of a network of oscillators described by

 ˙θj=ωj+λN∑n=1Ajnsin(θn−θj). (1)

Due to the symmetric coupling, the system frequency is constant. In literature, natural frequencies are commonly shifted to fulfill without loss of generality, entering a reference frame that co-rotates with the system frequency. Moreover, one usually assumes that the network is connected, i.e., that along network links, each oscillator can be reached from any other.

If all natural frequencies are finite and the coupling strength surpasses a critical threshold , all oscillator frequencies eventually synchronize [7]; [8], which in the co-rotating reference frame is described by stable steady states of Eq. 1 where all oscillator phases are “locked”. Each steady state of Eq. 1 then belongs to a continuum of fixed points that differ only by an arbitrary uniform rotation of all oscillator phases. The onset and form of this frequency synchronization is characterized by the critical threshold and the values of locked phases, respectively, both of which depend on (the network topology) and the sequence of natural frequencies. While there are other flavors of synchronization that have been extensively explored, frequency synchronization is the regime of interest in many realistic settings where sparsely coupled oscillators have non-identical natural frequencies [3]; [8].

As an example, frequency synchronization is paramount to stable operations of transmission grids, i.e., power grids that predominantly operate on high voltages and alternating currents (AC) [3]; [9]; [10]. Knowing and predicting the critical threshold and phase-locked solutions enables a cost-effective grid design and modification of transmission grids [11]. Moreover, vulnerable transmission lines can be identified along which the power flow may exceed critical limits [12]. While voltage instabilities have been taken into account [13], assuming constant voltages suffices for all relevant purposes [10]. The AC power flow in transmission grids can be modeled with second-order Kuramoto-type equations [9], with steady states and their linear stability still captured by the original Eq. 1 [14]; [15]. To this end, it is legitimate to focus on a bipolar distribution of natural frequencies, partitioning the oscillator ensemble into power-injecting generators with uniform positive natural frequency and power-drawing consumers with uniform negative natural frequency. In this setting, consumers are commonly taken to represent substations that delegate power to lower-voltage distribution grids [9]; [10].

For many dynamical processes, the influence of their topological background has been thoroughly investigated (see [16] for a review); yet in the case of frequency synchronization described by the Kuramoto model, results on the effect of network structure are few (see reviews [2]; [5]; [16]). In [17], the exact critical threshold is implicitly computed for fully connected graphs and arbitrary natural frequency distributions. For the same setting, [18] delivers approximate analytic expressions for locked phases, yet does not capture those regimes just above critical thresholds for which locked phase differences exceed . For the same restricted class of steady states, a synchronization condition for arbitrary graph topologies is presented in [14] that is exact for several extremal topologies. In [19], the critical threshold is numerically determined for bipolar natural frequency distributions and various random network topologies. In [20], two highly regular ring-lattice configurations with symmetric bipolar frequency distributions are considered, obtaining lower bounds for the critical threshold through recording and evaluating mesoscopic network structures. These findings are extended to yield an approximative semi-analytical description for more heterogeneous network topologies, yet the stability of phase-locked solutions is not considered. In the collective-coordinate approach of [21], an educated guess of the functional form of locked phases is used to reduce the complexity of Eq. 1, obtaining analytically tractable evolution equations of the system (see [11] for further applications). This approach can be adapted to accommodate different coupling structures, yet requires the functional form of locked phases to be known a priori. To the best of our knowledge, there exists no analytic approach to date that, for power grid-like oscillator networks, predicts critical thresholds and locked phases other than for quite restrictive assumptions on graph topology, frequency mixing or the form of locked phases.

To fill this gap, we present in Sec. II a mean-field ansatz inspired by [22] and resembling the active-neighborhood approach previously used in the context of disease spreading [23]. In Sec. III, we use the ansatz to obtain a parametrized description of the effect of connection structure on frequency synchronization in oscillator networks - also for regimes just above the critical threshold. In Sec. IV, we present a simple self-consistent calculation of the respective parameter. We focus on the simple case of asymmetric bipolar distributions of randomly assigned natural frequencies, a more realistic setting than the symmetric and ordered setup in [20] to investigate stable operations of transmission grids. To characterize connection structure, we consider two simple local measures: an oscillator’s neighborhood size (its degree) and neighborhood composition (with respect to natural frequencies). This enables us to assess analytically the influence of an oscillator’s neighborhood on its locked phase, while allowing us to depart from the regular-lattice assumption in [20]. We assess results in Sec. V for regular random graphs where links are randomly distributed under the constraint that all oscillators have the same degree, and in Sec. VI for random graphs without aforementioned constraint. Our model captures phase-locking on configuration-model topologies [16] - simple network models that impose a given degree distribution on the network, i.e., a probability distribution of oscillator degrees, while otherwise featuring random connections.

## Ii Setting up mean-field equations

To model synchronization in transmission grids, we assume an integer number of generators with and , and an integer number of consumers with and . Hence the system frequency in Eq. 1 is zero, and steady states are identified with frequency synchronization in the system. Without loss of generality, we can constrain ; this accounts for the observation that real-world transmission grids usually contain more consumers than generators. For numerical investigations, we will mainly choose , as this is a realistic choice for the fraction of generators in transmission grids [10]. Moreover, all oscillator phases will be confined to in the following. In the spirit of [22], Eq. 1 can then be rewritten as

 ˙θj=ωj+λ[rjGsin(ΨjG−θj)+rjCsin(ΨjC−θj)] (2)

with and , multiplying both sides with in each case before equating the imaginary parts [6]. Hence the sum in Eq. 1 is split into two simple terms quantifying the coupling of oscillator to the group of adjacent generators and consumers, respectively. This comes at the cost of introducing four new variables per oscillator - one coupling amplitude and one neighborhood phase for each neighborhood type . Thus, in order to proceed analytically, we make the following simplifying assumptions:

(a1) We assume that the neighborhood phases and of any oscillator only depend on ’s natural frequency. This presupposes that the pull of ’s phase on neighboring oscillators’ phases is either (1) negligible, i.e., for sufficiently large neighborhoods or (2) dominated by , implying that is largely determined by . Hence for generators, we introduce , the mean phase of neighboring generators and , the mean phase of neighboring consumers [Figs. 1(a)-(b)]. In the same spirit, for consumers, we introduce as the mean phase of their neighborhood of generators and as the mean phase of their neighborhood of consumers.

(a2) For generator , we approximate the coupling amplitudes and as and , where and are the number of generator and the total number of neighbors of , respectively. This approximation supposes that there is a global mean field that faithfully reflects oscillator dynamics, and that an oscillator’s coupling to it is proportional to the oscillator’s local connectivity [22]. For well-connected, homogeneous graphs, this is the case [5]; [22] [Figs. 1(c)-(d)]. The proportionality factors and can be understood as the mean phase coherences of the generator and consumer portion of generator neighborhoods, respectively. Analogously, and can be defined, describing mean phase coherences in the two neighborhood types of consumers.

To compute for instance and , first consider that with (a1) and (a2), . Hence not only estimates the phase coherence of generator neighbors of a single given generator, but also measures the coherence of generator neighborhood phases across the whole generator ensemble. It therefore assesses the validity of assumption (a1), and does so the more accurately the better (a1)-(a2) are fulfilled. As without loss of generality, it follows that a value of close to indicates the goodness of (a1)-(a2), as exploited further below in assumption (a3).

With , it follows that . Similar expressions hold for , and , so that the coupling in Eq. 2 of an oscillator to the mean field just depends on , and . Consequently, all oscillators with the same natural frequency and neighborhood can be considered equivalent and described by a phase with [cf. Figs. 3(c)-(f) and Figs. 4(c)-(d)]. This is reminiscent of the active-neighborhood approach presented in [23] in the context of disease spreading, which captures configuration-model topologies characterized by the degree distribution . As frequency synchronization in a network presupposes the network to be connected, we set in the following. Assuming random mixing of natural frequencies, generators and consumers follow the same binomial distribution in and the same distribution in . Hence, with and , one obtains

 rYeiΨY = ∞∑k=1k∑x=0wg(k,x)gkeiθY(k,x) (3) ^rYei^ΨY = ∞∑k=1k∑x=0w1−g(k,k−x)(1−g)keiθY(k,x). (4)

In Eqs. 3-4, the time-dependent mean-field - initially defined through averaging over , i.e., over oscillator neighborhoods - is expressed through weighted averages over , i.e., over central oscillators. Note that in these weighted averages, oscillator classes contribute proportionally to their abundance and to the size of their relevant neighborhood type. This is in contrast to the classical definition of phase coherences that does not include neighborhood sizes in the weighting [6]. In the following, one can fix the value of each mean neighborhood phase coherence through a third assumption:

(a3) We set and , . As laid out above, this is consequential if one is convinced of the goodness of (a1)-(a2). More support for this assumption is given by the fact that the phase coherence of oscillators in a shared neighborhood is generally larger than this of oscillators of same type, but randomly picked in the network, since common neighbors are at most two links apart. As an increased intra-neighborhood phase coherence yields an increased inter-neighborhood phase coherence as quantified by and , (a3) follows [Figs. 1(c)-(d)].

With (a1)-(a3), Eq. 2 can be written as

 ˙θG(k,x) = 1+λxsin[ΨG−θG(k,x)]+λ(k−x)sin[ΨC−θG(k,x)] (5) ˙θC(k,x) = −g(1−g)−1+λxsin[^ΨG−θC(k,x)]+λ(k−x)sin[^ΨC−θC(k,x)]. (6)

Note that all information about network topology and frequency mixing is contained in the mean field of Eqs. 3-4, while the evolution Eqs. 5-6 for oscillator phases are unaffected otherwise. Multiplying both sides of Eq. 5 [of Eq. 6] with [with ], summing over and and adding both equations, one obtains with Eqs. 3-4

 ωS=−λkg(1−g)(^rG−rC)sin(^ΨG−ΨC) (7)

for as the system frequency. Note that in Eq. 7, and are not set to , but computed on the fly from Eqs. 3-4. Hence any inaccuracy in assumption (a3) will translate into a nonzero system frequency. Concurrently, it is easy to show along the same lines that Eq. 2 with just assumptions (a1)-(a2) (i.e., allowing for variable neighborhood phase coherences) has zero system frequency at any time . From the approximative character of (a3) follows that if the numerical integration of Eqs. 3-6 for sufficiently large does yield a synchronous regime (i.e., locked phase differences), it is generally not in steady state - despite the initial rescaling of natural frequencies to enter a co-rotating reference frame. Instead, all phases then synchronously rotate with the system frequency . If is subtracted from both sides of Eqs. 5-6, synchronization is reflected by a steady state of oscillators with natural frequencies shifted by and all phases at time rotated by .

From Eq. 7, is the smaller the better assumption (a3) is fulfilled. It vanishes for and , the latter due to the resulting symmetry in Eqs. 3-6 which yields . Yet also decreases with increasing and : both yields decreasing and a higher accuracy of the mean field, i.e., also decreasing . Consequently, in phase-locked regimes is largest for slightly supercritical and small , usually attaining values of less than one percent of generator and consumer frequencies. On the one hand, this small value justifies equating and in assumption (a3). On the other hand, it motivates a final follow-up assumption about the system:

(a4) One can neglect the small correction of the natural frequencies that arises when entering the co-rotating frame of Eqs. 3-6. This is achieved by fixing one arbitrarily chosen phase, yielding steady states whose form and linear stability is similar to those of synchronized solutions of the unaltered Eqs. 3-6.

In the following sections, the steady-state solutions of Eqs. 3-6 with assumption (a4) will be investigated and compared to the full system in Eq. 1, with stationary variables marked with a ‘*’-superscript. In all cases, numerical integration of Eq. 1 towards the phase-locked state is performed with the Heun scheme and step size until on connected networks, where sufficiently large network sizes are used to minimize correlations present in finite-size realizations of configuration models. To accurately model finite-size networks of size , all oscillator classes with an expected oscillator number smaller than should be disregarded in Eqs. 3-4, i.e., all classes for which in generator dynamics and for which in consumer dynamics. As , this necessitates a cutoff at at which for all , and at at which for all . Thus in the following, all sums over remain finite for finite network sizes, and is rescaled to fulfill .

The validity of assumptions (a1)-(a3) is then illustrated in Figs. 1(a)-(d) for generator neighborhoods in phase-locked states of regular random graphs. We choose this topology to specifically analyze how the mean field copes with heterogeneity in the neighborhood composition of oscillators. As expected, (a1)-(a3) become more accurate for larger degrees and coupling strengths. This is corroborated both by increasingly accurate averages and by decreasing variances for each oscillator class. Note that (a1)-(a3) already hold slightly above the respective critical threshold (red squares and blue diamonds), but also in Figs. 1(a)-(b) the slight dependence of mean neighborhood phases on the neighborhood composition. The latter implies a dependence of an oscillator’s phase on its neighborhood type that is investigated further below [cf. Figs. 3(c)-(f) and Figs. 4(a)-(b)]. Similar observations hold for consumer neighborhoods, thus additionally confirming (a4) via Eq. 7 (not shown). Furthermore, numerical integration confirms very similar critical thresholds and values of locked phases differences for Eqs. 3-6 with and without assumption (a4) (not shown). Thus with (a4), one can approximate synchronized solutions of the mean-field Eqs. 3-6 with its steady-state solutions for slightly altered natural frequencies. This makes the system amenable to analytic treatment, as laid out in the following.

## Iii Parametrized locked generator phases

As a consequence of assumptions (a1)-(a3), the Kuramoto model in Eq. 1 is well-approximated by two sets of neighborhood-class Eqs. 3-6 - one set for each oscillator type. These two sets can be disentangled through assumption (a4) and an appropriate choice of coordinates: with , and , Eqs. 3 and 5 become

 ˙θΨ(k,x) = 1−λ{xsin[θΨ(k,x)−Ψ]+(k−x)sin[θΨ(k,x)]} (8) Ψ = arg[kM∑k=kmk∑x=0wg(k,x)eiθΨ(k,x)]. (9)

Consequently, generator phases are just coupled to a single mean-field variable . This system co-evolves self-containedly and will be analyzed in the following; consumer dynamics can be decoupled and dealt with analogously (see further below).

A necessary condition for global phase-locking is that all generator classes phase-lock, which is surely the case if fixed points of Eqs. 8-9 are linearly stable for all . The complexity of this task can be reduced by first considering a constant parameter in Eq. 8 that is self-consistently computed via Eq. 9. This decouples the dynamics of each generator class not just from consumers, but also from other generator classes. The above phase-locking condition then translates into searching for values of which (1) yield a linearly stable fixed point in Eq. 8 for all (2) fulfill Eq. 9 and (3) yield a linearly stable fixed point in co-evolving Eqs. 8-9. The solution to the generator phase-locking problem can thus be divided into three parts: describing the steady state of Eq. 8 parameterized by a fixed (Appendix A) while characterizing some of its general properties (Appendix C), self-consistently computing with Eq. 9 (Appendix B), discarding unstable solutions (Appendices A-B) along the way .

Setting constant without specifying its value, vital properties of locked generator phases can already be analytically derived. Appendix A identifies and discards linearly unstable phase-locked states in Eq. 8, and calculates the closed-form expressions for linearly stable locked generator phases as

 sin[θ∗Ψ(k,x)] = λ−1[k−(1−cosΨ)x]+xsinΨ√k2−2x(k−x)(1−cosΨ)−λ−2k2−2x(k−x)(1−cosΨ) (10) cos[θ∗Ψ(k,x)] = sin[θ∗Ψ(k,x)][k−x(1−cosΨ)]−λ−1xsinΨ, (11)

with . Appendix B simplifies the self-consistent computation of in Eq. 9 to

 0=kM∑k=kmk∑x=0wg(k,x)sin[θ∗−Ψ(k,k−x)]≡FG(Ψ) (12)

with

 Ψ∈⎧⎪⎨⎪⎩[0,arccos(2[λkm]−2−1)][0,arccos(max{−1,2(λkm)−2−1−21−(λkm)−2k2m−1})] (13)

in case the smallest realized generator degree in the network is even- or odd-valued, respectively, and for odd-valued being reasonably asserted. Appendix C identifies universal monotonic behavior of locked phases; it shows that on the above interval,

 θ∗Ψ(k,x+1)≥θ∗Ψ(k,x) (14)

and

 θ∗Ψ(k′,k′x/k)≤θ∗Ψ(k,x) (15)

for and integer as well as

 θ∗Ψ(k,x)∈[arcsin(λk)−1,Ψ+arcsin(λk)−1] (16)

for stable locked generator phases, regardless of chosen degree distribution . Lastly, the ensemble phase spread - the maximum difference between two locked generator phases - is shown to be exactly , so that its upper bound is trivially given by Eq. 13. This upper bound becomes the tighter the closer the coupling strength is to the critical threshold. This is because computed stable steady-state in Eq. 12 are the closer to the upper bound of the search interval in Eq. 13 the smaller the coupling strength [Figs. 2(b)-(d)].

## Iv Determining parameter and locked consumer phases

To actually determine stable steady-state , decoupled generator dynamics with a co-evolving mean-field are captured in two steps: first, the parameter is self-consistently computed through Eq. 12 on the interval given by Eq. 13. Second, a linear stability analysis of co-evolving Eqs. 8-9 is performed at the fixed points given by self-consistent and the locked generator phases they parametrize in Eqs. 10-11 (cf. Appendix B). In Figs. 2(a)-(d), self-consistent values of - the roots of (red solid line) - are computed for a regular random graph. For sufficiently small , no self-consistent solution for exists [Fig. 2(a)] - there are no phase-locked generators. At , two solutions arise [Figs. 2(b)-(c)]. The larger-valued solution is computed to be linearly unstable and discarded. It disappears for sufficiently large to leave only one (small-angle) solution for that is always linearly stable [Fig. 2(d)]. Consequently, is the critical coupling strength in the decoupled generator ensemble, and generator phases can be stably locked in only one configuration. With in the unique stably locked regime at hand, - the steady-state mean phase of generator neighbors of a consumer - is obtained by inserting Eqs. 10-11 into Eq. 4 [red dashed line in Figs. 2(b)-(d)]. Consequently, analyzing the decoupled generator dynamics in co-evolving Eqs. 8-9 delivers three things: first the generators’ critical threshold , second the generators’ linearly stable locked phases , third the steady-state mean neighborhood phases and . All phases are given in coordinates in which generator natural frequencies are positive and , the latter as imposed in Eqs. 3-6 to decouple the two oscillator ensembles and obtain generator Eqs. 8-9.

Consumer dynamics can be described in the same self-contained manner by setting and in Eqs. 4 and 6. This is possible because the asymmetry of the two oscillator ensembles only lies in and in the resulting differing absolute value of their natural frequencies. Thus substituting in Eqs. 10-11 and 13-16 yields the respective quantities for the consumer ensemble if is now considered the number of consumer neighbors. With this change of parameters and indices, all considerations in Appendices A-C also apply to the decoupled consumer ensemble. The change in neighborhood indexing moreover leads to

 FC(Ψ)≡kM∑k=kmk∑x=0w1−g(k,x)sin[θ∗−Ψ(k,k−x)] (17)

(with in , cf. Eq. 10) for the function whose roots yield the self-consistent mean-field for the consumer ensemble. Numerically, one finds for the consumer ensemble in the regular random graph that for all considered and [see blue solid line in Figs. 2(a)-(d)], so that the critical threshold of the full system is given by . This can be qualitatively understood by considering that the lower bound on the critical threshold for generators is larger than that for consumers by a factor of (cf. Eq. 22 in Appendix A). Upon computation of the self-consistent stable - the steady-state mean phase of consumer neighbors of consumers - all consumer phases are given in a reference frame in which their mean generator-neighborhood phase is zero and their natural frequency is of positive value. For the stable steady state, the locked phases and [blue dashed line in Figs. 2(b)-(d)] are then computed analogously to the generator ensemble.

To relate locked generator phases [] and locked consumer phases [], one should express them in common coordinates. This can be done by (1) substituting in (2) changing the sign of all consumer phases (3) additionally rotating all consumer phases by an angle so that and . Step (1) expresses consumer classes based on generator neighbors, while step (2) accounts for the differing signs of the two natural frequencies. Note also that steps (1) and (2) leave Eq. 14 valid for the locked consumer phases in generator coordinates, whereas Eq. 15 becomes

 θ∗C(k′,k′x/k)≥θ∗C(k,x) (18)

for and integer . Finally, step (3) relates the two sets of locked phases by demanding that computing the same phase differently should yield (approximately) the same result. As a result, the rotation angle is over-determined as and , and the extend to which both equations yield similar gives a measure of the quality of approximations (a1)-(a4) that led to Eqs. 8-9. Figures 2(b)-(d) exemplarily show a good match already for weak supercritical coupling strengths in the case of the regular random graph, complementing Figs. 1(a)-(d) in supporting the validity of our mean-field assumptions. With the absence of bistability in both decoupled oscillator ensembles, one can also exclude for that type of regular random graphs the existence of bistable phase-locked regimes in the full mean-field description of coupled generator and consumer dynamics.

## V Comparing results for regular random graphs

With the mean-field framework fully laid out, its critical thresholds and locked phases can finally be computed through Eqs. 10-12 and their consumer counterparts. This is in turn compared with numerical integration of the full system in Eq. 1. To illustrate the role of neighborhood composition in phase-locking, we first choose regular random graphs [sparse graphs in Figs. 3(a), (c), (e), denser graphs in Figs. 3(b), (d), (f)]. As in Figs. 1(a)-(d), we immediately notice decreasing intra-class variances for increasing graph connectivity and coupling strength. This confirms that our nearest-neighbor approach works best for well-connected graphs where longer-range correlations cause little variability of locked phases within the same oscillator class.

For the critical thresholds in Figs. 3(a)-(b), one observes that the closer is to , the larger is the critical threshold. This is plausible, as then the difference in natural frequencies to bridge for synchronization increases, while the phase pull of the weaker generator ensemble also increases due to the ensemble’s increasing size. Our framework (i) systematically underestimates ; deviations grow with (ii) decreasing and (iii) increasing . The reason for (i) is that our ansatz only distinguishes oscillators based on their natural frequency and neighborhood type, yielding a mean-field description for the phase-locking of each oscillator class. Yet the onset of frequency synchronization in real networks presupposes that also outlier oscillators phase-lock - oscillators whose network embedding (i.e., connections beyond their neighborhood) is particularly detrimental to their synchronization with respect to the mean-field description of their class. These outliers determine the onset of synchronization in real networks, retarding it with respect to our mean-field approach. This retardation is particularly pronounced for sparse connectivity, as there fluctuations in the topology and frequency composition of higher-order neighborhoods are more manifest in central oscillators’ dynamics, thus explaining (ii). Finally, higher-order frequency-composition fluctuations with retarding effect are more pronounced the closer is to , accounting for (iii). This is because, combinatorially speaking, there are more higher-order neighborhood configurations possible for more similar ensemble sizes. A higher number of configurations also means a higher chance of fluctuations with retarding effect on phase-locking. For the same reason, the retarding effect of outliers is more pronounced in topologies with broader degree distributions like the random graphs dealt with in Sec. VI.

For increasing supercritical coupling strengths, the importance of said fluctuations diminishes [compare Fig. 3(c) with (e) as well as Fig. 3(d) with (f)]. This becomes apparent first through the increasingly accurate prediction of locked phases by our framework (i.e., by considering Eqs. 10-12, 17 and 3-4 in that order), second in the full system through the decreasing variances of locked phases within one oscillator class. Moreover, plotting in Figs. 3(a)-(f) the output of Eq. 2 with only assumptions (a1)-(a2) yields very similar results to our final approach that also includes (a3)-(a4) [dashed lines in Figs. 3(a)-(b), complete overlap with solid lines in Figs. 3(c)-(f)]. This underlines the quality of (a3)-(a4) and indicates that (a1)-(a2) are the most significant source of deviation from the full system given by Eq. 1. As predicted in Eq. 14, locked phases increase monotonously with , and, for fixed natural frequency, spread no farther apart than predicted by Eqs. 13 and 16 and their consumer counterparts. We note that for coupling strengths just above the critical threshold, the ensemble and thus also the global phase spread can exceed [Figs. 3(c)-(d)]. This is a phase-locked regime commonly not considered in literature where instead, authors build on a result from [14] that shows the linear stability of all steady states with global phase spread smaller than . Our model is also in line with this result: observed unstable self-consistent solutions of Eq. 12 for the regular random graph all possess global phase spreads larger than .

## Vi Application to other connection structures

In highly regular random graphs with the same neighborhood composition for all oscillators, (an integer) is the only relevant oscillator class. Hence with the Dirac function, so that Eqs. 3-4 imply at all times for oscillator type . Therefore the system reduces to an effective two-oscillator problem. For both decoupled generator and consumer dynamics, Eqs. 10 and 12 (and their consumer counterparts) yield for the respective free parameter. In both cases, this parameter is trivially self-consistent as then . Hence and also for locked neighborhood phases, so that locked consumer phases can be unambiguously expressed by generator coordinates. In these coordinates, the difference between locked generator and consumer phases is then the relevant system variable. Its linear stability analysis in Appendix B then reduces to determining when . Considering that and , this reveals a unique stable state at above the critical threshold .

The same results are obtained by a more naive mean-field approximation of Eq. 1. There, oscillators of same natural frequency can be considered equivalent through setting for and for . Splitting the coupling term into two contributions from interactions with generators and with consumers, as well as averaging Eq. 1 over the generator and consumer ensemble yields for the phase difference , with steady states and stability as above. Hence our framework is consistent with a more naive mean-field approximation, which in our setup moreover coincides with the collective-coordinate approach [21] for regular random graphs. In Figs. 3(a)-(f), the predictions of the naive mean-field model are plotted with dotted lines, showing that the inclusion of neighborhood heterogeneity in our more advanced approach yields a much better approximation of Eq. 1.

The application of the mean-field approach to connected random (Erdős-Rényi) graphs with is illustrated in Figs. 4(a)-(f). For the chosen mean degree and network size , the minimum considered degree is and maximum degree (cf. Sec. II), resulting in a cutoff for larger degrees [extended white area in Figs. 4(e)-(f), cf. Sec. II]. With and Eq. 23, a necessary condition for global phase-locking is . Choosing , we find the full system to be in a weakly supercritical regime in which we test the validity of our mean-field approach.

We observe in Figs. 4(a)-(b) that the monotony of oscillator phases in the full system of Eq. 1 is correctly predicted by Eqs. 14, 15 and 18. Furthermore, locked oscillator phases converge to the reference phase for large degrees. This is because the larger an oscillator’s neighborhood is, the stronger is its coupling to the rest of the network and the mean-field variable describing it. That also explains the convergence of locked consumer phases towards for increasing degrees [Fig. 4(b)], as this is the mean phase consumers couple to in their ensemble description. As in the case of regular random graphs, the overall phase spread exceeds . For the full system, the variances of locked phases within an oscillator class is very low for both oscillator ensembles [Figs. 4(c)-(d)], also for classes describing oscillators with lower degrees. This validates our mean-field approach of only considering oscillator neighborhoods when modeling phase-locking in configuration models. We normalize by the respective ensemble spread to obtain the relative magnitude of fluctuations. Lastly, Figs. 4(e)-(f) illustrate the differences of averaged locked phases in the full system and predicted locked phases, again normalized by the respective ensemble spread. They show that our framework predicts well the locked phases in the full system already for weak supercritical coupling strengths. As in Figs. 3(c)-(f) for regular random graphs, we find a very good match for oscillator classes ) for which , while predicted locked phases deviate the more the larger are deviations of from that mean number of generator neighbors.

## Vii Summary and conclusions

We present a mean-field approach [Eqs. 3-6 with assumption (a4)] to analytically assess how the phase-locking of randomly coupled oscillators is affected by their local connection structure. To that end, we consider asymmetric bipolar distributions of natural frequencies, modeling the most common asymptotic regime in electrical transmission grids. Generator and consumer ensembles can be decoupled, each of which is analyzed separately and analogously (Eqs. 8-9). For each ensemble, unstable phase-locked regimes are detected and discarded, and stably locked phases are given closed-form parametrized expressions (Eqs. 10-11). The respective functional form of locked phases is independent of imposed degree distributions and frequency mixing. Such global information about the network is instead contained in the free parameter that acts as a mean field to which the considered oscillator ensemble couples (Eq. 9). For each ensemble, the free parameter can be self-consistently computed through a simple expression (Eqs. 12 and  17, respectively). The search space can be analytically constrained (Eqs. 13), yielding simple statements about the locked phases’ bounds (Eq. 16) and monotony (Eqs. 14, 15 and 18) already without knowing the parameter’s numerical value. As a result, our framework predicts the same functional form and monotony of locked phases for configuration-model networks of any degree distribution, provided that the underlying mean-field assumptions are fulfilled reasonably well.

Calculating the self-consistent free parameter, we find that above an ensemble-specific critical threshold, there exists a linearly stable steady state in each decoupled ensemble dynamics. For coupled ensemble dynamics, i.e., the full mean-field model, this automatically yields the overall critical threshold. Furthermore, for examined parameters and topologies, the framework rules out coexisting linearly stable phase-locked regimes that conform to mean-field assumptions (a1)-(a4). In order to couple ensemble dynamics, a simple coordinate transformation of locked consumer phases is needed given by an overdetermined rotation angle. The discrepancy of its calculated values is a model-intrinsic measure for the accuracy of used mean-field assumptions, as is (the absolute value of) the system frequency of the mean-field equations (Eq. 7). Both measures indicate the validity of our mean-field approximations already for sparse and weak coupling. Additionally, we find that our approach is consistent with a simpler mean-field ansatz, lending further support to the used mean field.

Naturally, the framework’s accuracy must ultimately be assessed through comparison with the full system: for two configuration-model topologies and already for weak coupling, we observe that closed-form locked phases match well averaged phases in the full system. We moreover find that for coupling strengths just above the critical threshold, the largest difference of locked phases exceeds in both the full system and its mean-field approximation - a regime commonly not considered in literature [3], but captured by our approach. For densely connected regular random graphs, critical thresholds are well-predicted; more accurately than through an alternative approach in [21]. Their systematic underestimation there as well as in sparser and more heterogeneous topologies is tied to the importance of outlier oscillators that do not obey mean-field assumptions.

In future work, a more systematic investigation with our framework could reveal topologies with coexisting phase-locked regimes. While such multistability has been detected for single network realizations [24], our framework can potentially detect multistable regimes of a whole class of networks given by the respective configuration model. Furthermore, by adjusting the binomial weights in Eqs. 3-4, our assumption of random frequency mixing could be relaxed to account for frequency correlations as in [25]. Finally, by yielding the functional form of locked phases a posteriori, our proposed approach could moreover tie dimension-reduction approaches such as in [21]; [26]; [27] to intuitive mean-field assumptions.

###### Acknowledgements.
The authors would like to thank F. Daviaud, F. Suard, M. Velay and M. Vinyals for useful discussions. Funding by CEA under the NTE NESTOR grant is gratefully acknowledged. Additionally, S.W. acknowledges funding of project 330 in the Enhanced Eurotalents program.

## Appendix A Parameterized locked generator phases

Here, we want to investigate the existence and form of stable phase-locked solutions of Eq. 8 for a given degree . To this end, we set constant in the following and rescale time by a factor , so that the decoupled generator dynamics of Eq. 8 can be rewritten as

 ˙θΨ,k(z)=(λk)−1−fΨ(θ,z) (19)

with , , , integer and . The contour lines defined by are called solution branches in the following, because if they encompass , then yields a steady state in Eq. 19 and generator class is phase-locked in Eq. 8 for fixed (stably or unstably). To understand how solution branches arise depending on and for given , consider moreover the maxima of curves for each given by

 ^fΨ(z)≡√1−2z(1−z)(1−cosΨ). (20)

Clearly a solution branch will encompass as soon as

 ^fΨ(z)≥(λk)−1. (21)

A short inspection of confirms that if is a steady state, so are and all for which , with the necessary coupling strength increasing the closer is to . Hence if oscillator class phase-locks, so do all classes with and , with the class of the most balanced numbers of generator and consumer neighbors phase-locking last.

Obviously always holds, so that for , no solution branch exists for any . In this case, there is no steady state in Eq. 19 for any or , so that no oscillator class in Eq. 8 phase-locks. For larger , two joined solution branches and appear, related as and confined to all fulfilling . Hence only generators with a sufficiently small or large number of generator neighbors phase-lock. Finally, for all coupling strengths larger than

 λ∗Ψ,k=[k|cos(Ψ/2)|]−1, (22)

both solution branches are separated, and each lives on the entire interval . As a consequence, all generator classes for a given phase-lock in Eq. 8, in particular class in graphs with even-valued degree . In co-evolving Eqs. 8-9, (computed with steady-state ) obviously is a lower bound for the critical threshold if is even-valued.

For an odd-valued degree , this lower bound is approximate: there, only the two generator classes at and must lie on the solution branches for all classes in Eq. 8 to phase-lock, which is the case as soon as

 λ≥[(λ∗Ψ,k)−2+sin2(Ψ/2)]−1/2 (23)

(cf. Eqs. 20-22). This lower bound is always smaller than , and the better approximated by it the larger and the smaller are. Their largest relevant difference is obtained by maximizing on and minimizing : in Eq. 12 and Appendix B, we observe to be the approximately largest self-consistent stable value for various degree distributions and generator abundancies . Considering that is the minimum odd-valued degree, one calculates versus the exact lower bound in Eq. 23. Choosing the next-largest odd-valued degree , we already obtain versus the exact lower bound . Conversely, for the actual stable phase-locked regimes in co-evolving Eqs. 8-9, we observe in Eq. 12 and Appendix B, again for various and . Hence in the following, we will assume for odd-valued degrees that for stable steady-state in co-evolving Eqs. 8-9. This is a weaker statement than the strict proof for even-valued that for any in phase-locked states of Eq. 8.

Identifying [] with the solution branch for which [] for , phase-locked states lying on are unstable and can thus be discarded: firstly, observe that any linearly unstable steady state in Eq. 8 with fixed cannot be a stable phase-locked state in the system of Eqs. 8-9 where co-evolves. In turn, for oscillator classes with even-valued (odd-valued) degree , Eq. 22 (Eq. 23) is a necessary condition for phase-locking in co-evolving Eqs. 8-9. Secondly, the Jacobian for Eq. 19 with fixed and is and thus changes sign at the maximum of . As the line of these maxima separates the two solution branches (Eq. 20) and furthermore holds, it follows that the Jacobian is always positive for fixed points on [and negative for all fixed points on ]. Hence the fixed points of Eq. 19 are unstable and discarded, and the stable steady states are given by . They are of the form

 sin[¯θ∗Ψ,k(z)] = (λk)−1[1−(1−cosΨ)z]+zsinΨ√1−2z(1−z)(1−cosΨ)−(λk)−21−2z(1−z)(1−cosΨ) (24) cos[¯θ∗Ψ,k(z)] = sin[¯θ∗Ψ,k(z)][1−(1−cosΨ)z]−(λk)−1zsinΨ, (25)

with . From the behavior of discussed above, it furthermore immediately follows for the contour lines that

 ∂¯θ∗Ψ,k(z)/∂(λk)≤0. (26)

## Appendix B Self-consistent generator mean field

A parameterized phase-locked solution of Eq. 8 describes the same phase-locked state as , but in other angular coordinates, namely for and in Eq. 5 instead of and as before. This coordinate shift should not change the physics of generators in the system; it can however be used to simplify the computation of self-consistent in Eq. 9. As its left-hand side is just , Eq. 9 reads

 0=kM∑k=kmk∑x=0wg(k,x)sin[θ∗−Ψ(k,k−x)]≡FG(Ψ) (27)

in the shifted coordinates, to be fulfilled by all self-consistent stationary .

A generator class with degree curtails the search space for these through Eqs. 22-23: with necessary for any class with degree to phase-lock (Eqs. 19 -21), all self-consistent must fulfill

 cosΨ∈⎧⎪⎨⎪⎩[2(λk)−2−1,1][max{−1,2(λk)−2−1−21−(λk)−2k2−1},1] (28)

for even- and odd-valued degree , respectively. The search space of self-consistent is thus dictated by the smallest realized generator degrees.

To determine the linear stability of self-consistent phase-locked solutions of Eqs. 8-9 with co-evolving mean field , first consider that through Eq. 9 and ,

 ˙Ψ=(gkrG)−1kM∑k=kmk∑x=0wg(k,x)˙θΨ(k,x)cos[Ψ−θΨ(k,x)],

with as in Eq. 3. The latter reappears after having set already in assumption (a3). Methodologically, it is coherent to re-set , at the cost of rendering the following stability considerations inexact. Next, one should determine the entries of the system’s Jacobian, computed at the stationary state parametrized by . For , and , these are, at the fixed points,

 ∂˙θΨ(k,x)∂θΨ(k′,x′)∣∣θ∗ = −δkk′δxx′λ√k2−2x(k−x)(1−cosΨ)−λ−2 ∂˙θΨ(k,x)∂Ψ∣∣θ∗ = λxcos[Ψ−θ∗Ψ(k,x)] ∂˙Ψ∂θΨ(k,x)∣∣θ∗ = wg(k,x)gkrGcos[Ψ−θ∗Ψ(k,x)]∂˙θΨ(k,x)∂θΨ(k,x)∣∣θ∗ ∂˙Ψ∂Ψ∣∣cθ∗ = λgkrGkM∑k=kmk∑x=0xwg(k,x)cos2[Ψ−θ∗Ψ(k,x)]

with

 cos[Ψ−θ∗Ψ(k,x)]=[x+(k−x)cosΨ]√k2−2x(k−x)(1−cosΨ)−λ−2+λ−1(k−x)sinΨk2−2x(k−x)(1−cosΨ)

from Eqs. 10-11. Consequently, the Jacobian’s eigenvalues are also parameterized by ; inserting the latter’s self-consistent solutions computed in Eq. 12 reveals the linear stability of the steady states of co-evolving Eqs. 8-9. If all computed eigenvalues have negative real parts, the system is linearly stable. If at least one eigenvalue possesses a positive real part, the system is linearly unstable.

## Appendix C Properties of locked generator phases

First, we show for fixed even-valued degree that if , then