Graphop Mean-Field Limits and Synchronization for the Stochastic Kuramoto Model

Models of coupled oscillator networks play an important role in describing collective synchronization dynamics in biological and technological systems. The Kuramoto model describes oscillator's phase evolution and explains the transition from incoherent to coherent oscillations under simplifying assumptions including all-to-all coupling with uniform strength. Real world networks, however, often display heterogeneous connectivity and coupling weights that influence the critical threshold for this transition. We formulate a general mean field theory (Vlasov-Focker Planck equation) for stochastic Kuramoto-type phase oscillator models, valid for coupling graphs/networks with heterogeneous connectivity and coupling strengths, using graphop theory in the mean field limit. Considering symmetric odd-valued coupling functions, we mathematically prove an exact formula for the critical threshold for the incoherence-coherence transition. We numerically test the predicted threshold using large finite-size representations of the network model. For a large class of graph models, we find that the numerical tests agree very well with the predicted threshold obtained from mean field theory. However, the prediction is more difficult in practice for graph structures that are sufficiently sparse. Our findings open future research avenues toward a deeper understanding of mean-field theories for heterogeneous systems.


I. INTRODUCTION
The discovery of synchronization dates back to 1665 with Christiaan Huygens' observations of two synchronizing pendulum clocks 9 , and its mathematical modeling likely began with Norbert Wiener who was inspired by neuronal oscillations in the brain 10 .Wiener's formulation of the problem, a) Corresponding author: erik.martens@math.lth.sehowever, was too general to allow for any analytical progress; simplifying assumptions were necessary to render the problem mathematically tractable 11 , culminating in Yoshiki Kuramoto's paradigmatic model 4,5 .Kuramoto's original model describes the time evolution of the oscillator phases θ k = θ k (t), where k ∈ {1, . . .N} =: [N], the coupling interaction between oscillators is first order, the coupling is all-to-all with uniform strength C, and the intrinsic frequencies ω k are drawn unimodally from a distribution g centered in the origin.The level of synchronization in this transition is aptly captured using the order parameter r(t) = 1 N ∑ N j=1 exp (iθ j (t)) , which tends to 0 when oscillator phases are incoherent (disordered) for weak coupling or to 1 when oscillators lock their frequencies and phases clump together (we say that the phases are coherent / the oscillators synchronize).When frequencies are identical, ω k = ω j for all k, j ∈ [N], the so-called synchronization manifold defined by θ k (t) = θ j (t) for all k, j ∈ [N] exists and is attractive for C > 0; vice versa, when frequencies are nonidentical (or symmetry is broken due to some other mechanism, see below), the loss or gain of coherence plays out in a competition between the strength of the heterogeneity and coupling strength C. Thus, for a set distribution width of the intrinsic frequencies, Kuramoto's model exhibits a transition from incoherent to coherent oscillations as the coupling strength C surpasses a certain threshold value C .Kuramoto's initial heuristic analysis was based on a selfconsistency equation for the order parameter 12 , allowing one to predict the critical coupling strength associated with the incoherence-coherence transition.A more formal mathematical treatment, facilitating deeper insights would, however, require a mean field theory valid in the limit, N → ∞.Such a theory 13 describes the dynamics in terms of a density function in the oscillator phases, ρ = ρ(φ ,t), which evolves according to a transport equation (formally, a Vlasov-Fokker-Planck equation; see Eq. ( 8)).Such a description was used by Strogatz and Mirollo 14 to investigate the stability of the incoherent branch for C < C where ρ = 1/(2π) by studying the associated eigenvalue spectrum and to (re-)derive the critical coupling C = 2/(πg(0)), where g(0) denotes the maximum value of a unimodal frequency distribution g; this approach was further developed and applied to variants of the Kuramoto model 15 .Other studies focused on the stability analysis for the partially synchronized branch (C > C ) 16 .An exact low dimensional description in terms of the macroscopic dynamics (order parameter) allowing to express the evolution of the order parameter in terms of an ordinary differential equation became available later [17][18][19] .While Kuramoto's simplifying assumptions allowed for making significant progress in the mathematical understanding of the synchronization phenomenon, to understand real world oscillator dynamics, it is desirable to break these assumptions toward increasing complexity.There are a number of ways of doing this; here, we are concerned with how the incoherence-coherence transition is affected by the presence of (thermal) noise and, in particular, network heterogeneities, which play a major role in real systems [20][21][22][23][24][25][26][27] .Indeed, the ability of coupled oscillators to synchronize has been investigated under the influence of noise 12,28 , heterogeneous connectivity 29,30 , or heterogeneous coupling, such as non-local 31,32 , k-nearest-neighbor 33 , or random coupling strengths 34,35 and also on experiments [36][37][38][39] where oscillators are subject to real world influences.
Mean-field descriptions for N → ∞ are well established for various theoretical frameworks including coupled oscillator networks 13 .Our focus thus lies on mean-field limits valid for complex networks 8 , i.e., to generalize the Vlasov-Fokker-Planck (VFPE) equation (see Eq. ( 8)) so that it is capable of accurately describing the dynamics in complex networks characterized by heterogeneities in the connectivity or coupling strength.In particular, this includes cases where the adjacency matrix defining interactions between finitely many vertices is neither a full graph nor a highly symmetric structure, such as a lattice.In order to incorporate such structures, it is necessary to extend the description of (weighted) graph structures to the mean-field limit.This is possible via so-called graphons, which rely on concepts of the theory of limits of graph sequences 40,41 or even more generally utilizing the theory of graphops 42 .Intuitively, graph limit theory provides a way to arrange limits of discrete graphs as continuous objects.Graphons achieve this, mostly within the context of dense graphs, using a coupling kernel function that describes the connectivity in the limit.Graphops generalize graphons, also incorporating many intermediate and sparse density graph limits in addition.Graphops can be represented as operators or via an associated measure-theoretic representation; i.e., they are generalizing purely kernel-based operators to more general operators.Recent studies have used graph limit theories to pursue the goal of heterogeneous mean-field limits.Several mathematical approaches have been successful in providing rigorous proofs for VFPEs, where nonlocal integral terms appear to take into account the heterogeneous coupling structure 43,44 .Recently, a general theoretical framework based on graphops has been put forward (by some of the authors of this paper) that allows us to generalize mean-field limit VFPEs easily from particular cases (nonlocal coupling or standard all-to-all) to describe modern complex network structures [45][46][47][48] .
In the present paper, we extend previous work [45][46][47] to the stochastic case and formally derive a mean-field description based on graphop theory for the Kuramoto model with identical oscillators interacting via first order harmonic coupling function and non-uniform coupling strengths under the influence of (thermal) noise.We then derive rigorous results for the critical threshold for the incoherence-coherence transition (C ) by deriving a stability formula for the incoherent solution branch.A difficulty arises as it is unclear what demarcates the boundary of validity of mean-field PDEs for complex heterogeneous graphs; i.e., at some level of graph heterogeneity, it may be too difficult to accurately capture details of very sparse graph structures.As we cannot be sure under what circumstances our results correspond to the dynamics obtained for finite graphs (rigorous convergence results are still needed), we carry out detailed numerical simulations to test our results for various finite graph structures.
This article is structured as follows.In Sec.II, we introduce a formal derivation of the mean-field equations for N → ∞.In Sec.III, we derive the critical coupling strength C for the continuum limit, based on the graphop mean-field limit equation.In Sec.IV, we carry out numerical simulations to investigate how the incoherence-coherence transition point predicted by the mean-field theory carries over to finite graphs for a range of graph structures, including dense and sparse topologies.Finally, we discuss our results in Sec.V.

II. FORMAL DERIVATION OF THE MEAN-FIELD EQUATIONS
As discussed above, we are interested in mean-field models for stochastic Kuramoto(-type) models on networks 49 .The individual dynamics for the coupled identical oscillators is given by where k ∈ [N] := {1, ..., N} and θ k = θ k (t) ∈ T := R/(2πZ) denotes the phase of the kth oscillator, (A N k, j ) k, j∈[N] denotes a weighted and non-negative adjacency matrix of the network (i.e., a graph G with adjacency matrix A N k j ), D : T → R is a sufficiently regular coupling function (e.g., D = sin), C > 0 is the coupling strength, β > 0 is a diffusion constant controlling the noise level, and W (t) = (W 1 (t), . . .,W N (t)) is a vector of N independent Brownian motions so that Ẇk = Ẇk (t) is just a white noise forcing for each oscillator.The following derivation extends 45 to the stochastic case.To understand the formal derivation, let us consider the classical case of the Kuramoto model with all-to-all coupling with uniform strength, i.e., for A N k, j = 1 for all k, j ∈ [N] and D = sin.In other words, (1) now reads as Let us introduce the complex order parameter Multiplying this equation by e −iθ k and equating imaginary parts, we have which implies that From ( 5), the mean-field character of the problem is visible as the kth oscillator just feels the averaged input from all other oscillators so one can think of a single typical oscillator and aim to analyze its dynamics.Let ρ(t, θ ) dθ denote the fraction of oscillators with phase between θ and θ + dθ at time t; i.e., ρ is a probability density.Assuming a law of large numbers in the limit N → ∞, we formally get Now, using the same trick as above (i.e., multiplying both sides of the last equation by e −iu and taking imaginary parts), Eq. (5) becomes in the limit Finally, the continuity equation, also called the Vlasov-Fokker-Planck equation (VFPE), for the probability density ρ, respectively, for the law of the limiting process u, reads as In summary, ( 8) is a partial differential equation with a firstorder transport/advection-type term with a nonlocal convolution term involving the sine-nonlinearity mediating the coupling and with a second-order spatial diffusion term arising directly from the white noise forcing.Now, let us come back to Eq. (1).In this case, the next natural generalization step is to assume that the network (i.e., a graph G with adjacency matrix A N k j ) is sufficiently connected and does not have components, which are more connected than others; see also Ref. 49 .Moreover, we assume that there exists a local order parameter r k e iψ k , which is locally proportional to a single global order parameter re iψ weighted by the degree κ k for each node; i.e., we have By multiplying the local order parameter by e −iθ k and equating the imaginary parts in the last equation, we obtain Now, let ρ(t, θ , κ) dθ denote the probability for the fraction of oscillators having a phase between θ and θ + dθ and a degree κ at time t.Note carefully that we have added an additional variable κ to the density, which captures the (degree) heterogeneity of the network.If we assume that the network is uncorrelated and has degree distribution d(κ), one is tempted to assume that in the limit N → ∞, we have where κ is the average degree of a vertex in the graph and κd(κ) κ ρ(t, φ , κ) is the probability density for an edge having its end at a vertex of phase φ and degree κ at time t.Now, using the same trick as before, Eq. (10) becomes in the limit The continuity equation for the probability density ρ, respectively, the law of the limiting process u, reads as Thus, in comparison with the VFPE for the classical case of all-to-all coupling with uniform strength (8), we had to replace We can view this step as incorporating the structure of graph/network G appearing in the Vlasov equation via an operator, which acts on the density ρ.In fact, one can even hope to completely remove averaging over the variable κ that we used to capture the heterogeneity and just keep κ as a new variable in the density, which then yields a whole hierarchy of mean-field VFPEs, one for each degree.This set of ideas can then be thought even further and one can directly replace the adjacency matrix by a coupling kernel and there are numerous papers in this direction 43,44,46,47 .Yet, it seems best to think of generalizing VFPEs more abstractly 45 by viewing the underlying network influence as given by some linear operator A acting on the density so that a more abstract form of VFPEs would be given by where x is a suitable variable that tracks the heterogeneity of the network so that one effectively obtains a family of VF-PEs, and we have also replaced the sine-coupling again by a more general coupling function D. A typical choice of x found in the literature would be to take it as a variable in the unit interval x ∈ [0, 1] = Ω, where points in the interval represent node labels in the infinite network limit 43,44,46,47 .
Probably the most elegant abstract way to think of A is as a graph operator, or graphop, as introduced in Ref. 42 .A graphop is a bounded, self-adjoint, and positivity-preserving operator A : , where m is the reference measure on Ω; e.g., one can pick the Lebesgue measure.To a given graphop A always corresponds a family of finite measures (ν x ) x∈Ω , called fiber measures, via the formula Intuitively, we may view a graphop A just as a generalized adjacency matrix for a symmetric graph and for a given node x ∈ Ω, the fiber measure ν x is just the edge distribution for this node.Indeed, for the finite-dimensional case, we can just pick Ω = [N] and m as the uniform measure on Ω, so that functions f ∈ L ∞ (Ω; m) can be identified with vectors in R N and A f is just the usual matrix-vector multiplication.Yet, we stress that in the limit N → ∞, we need a space, such as Ω = [0, 1], with the Lebesgue measure.One may wonder, how far such an abstract construction for VFPEs involving graphops can work?It is clear that it works in simple cases, e.g., when the graph is all-to-all coupled as one can just drop the dependence on x.Also, if the graph is very dense and very regular with just two types of typical nodes, then one could take x as a binary variable and so on.Furthermore, it is understood that it works for dense graphs, where A can be represented by an integral operator with a sufficiently regular kernel, i.e., in the framework of so-called graphons.However, one does expect that there are growing networks as N → ∞ that are so sparse and/or so heterogeneous that eventually, mean-field calculations may fail.Proving a precise boundary location on the space of networks to determine, when VFPEs are helpful and when they fail, seems out of reach at this point.Here, we take a pragmatic approach and start from the formal VFPE (15a), carry out stability analysis of the main bifurcation/phase transition to synchronization in the Kuramoto model, and then numerically simulate the dynamics for different discretized (i.e., finite-dimensional, large N) classes of graphops A to check when the mean-field stability calculation is accurate.This is going to provide an indirect cross-check, whether a mean-field limit can work.

III. BIFURCATION/PHASE TRANSITION
In the following, we consider (15a), and we assume for simplicity that (H0) The coupling is non-trivial; i.e., C > 0.
(H2) A is a graphop with a bounded 2 → 2 norm; i.e., the following quantity exists and is finite: This implies that A can be uniquely extended to the Hilbert space L 2 (Ω, m) (see Remark 2.12 in Ref. 42, for instance).For simplicity, we use the same notation for this extension; i.e., we write A : L 2 (Ω, m) → L 2 (Ω, m).For the solution ρ(t, θ , x) of (15a), we define the jth Fourier coefficient as where i := √ −1.Note that we have effectively defined a family of Fourier coefficients that depends upon x, i.e., {(z j ) x } x∈Ω , but we shall always write just z j in the calculation below and later discuss the x-dependence.Applying the Fourier transform to (15a), exchanging integrals, and using integration-by-parts (in the second line), we have where D denotes the Fourier transform of D, and in the last line, we used that D(0) = 0, which follows from the fact that D is an odd, periodic function.Moreover, z − j = z j holds, which follows from the fact that ρ is real-valued.We can assume without loss of generality that j ∈ N to get the following system (i.e., the amplitude equation): The completely incoherent state ρ ∞ ≡ 1/(2π) of the oscillators corresponds to a uniform probability density over the circle, which translates into z 0 = 1/(2π) and z j = 0 for all j = 0, and the state is also assumed to be independent of x; i.e., we assume that all different types of nodes are uniformly distributed across the circle for ρ ∞ .Linearizing (18) around this incoherent state yields via a straightforward calculation the system where we use Z j to denote the Fourier coefficients of the linearized dynamical system, and we observe that the linearized system nicely decouples.The question then is how does the stability of the jth Fourier mode depend on the eigenvalues of graphop A? On the Hilbert space H := L 2 (Ω, m), for any j ∈ N, let us define the linearized operator Recall that the resolvent set ρ(A) of the operator A : H → H is defined to be the set H → H exists and is bounded}, where R λ (A) is called the resolvent operator of A and the spec- From this, we see that for all j for which D(− j) = 0, the condition that R λ (A) exists and is bounded is equivalent to the condition that R λ (T j C ) exists and is bounded.From this, we conclude that for all j ∈ Z for which D(− j) = 0, we have For all other j ∈ Z (that is, for all j for which D(− j) = 0) we see immediately that σ (T j C ) = − j 2 β 2π .Since A is bounded and self-adjoint, we have that σ (A) ⊂ R is a bounded set.Further note that since D is an odd function, we must have D( j) = i 2π 0 D(u) sin( ju)du ∈ iR.Finally, define where Z * := Z + ∪ {0}.The next theorem shows that C is a uniform parameter bound on the coupling strength independent of x, which means that smaller coupling leads to stability of incoherence, while above C , at least some classes of nodes synchronize at least partially.More precisely, we have Proof.Observe that for any λ ∈ σ (A) and j ∈ Z * , such that i D(− j)λ < 0, the corresponding element in the spectrum of ), is strictly negative for any C > 0; thus, it never crosses the imaginary axis.Thus, a crossing, for growing C, can occur only among those λ ∈ σ (A) and j ∈ Z * for which i D(− j)λ > 0. Among all such λ and j, the crossing occurs always at Observing that C is the minimum of all these transition points, it follows immediately that C is the smallest C > 0 for which there exists j ∈ N such that an element in the spectrum of T j C crosses the imaginary axis, namely, the element λ (C , j) ∈ σ (T j C ).As a first step, we want to carry out some specializations to examples and analytically consider some cases.

Example III.2. (Kuramoto model with first order interaction)
We wish to specialize the general formula for the critical threshold in Eq. ( 20), valid for Eq. ( 15) [the continuum limit version of Eq. ( 1)], to the continuum limit version of the Kuramoto model in Eq. ( 1) such that oscillators interact via a first-order coupling function, i.e., D(u) = sin u.Hence, and i D(−1) = 1 2 > 0.Then, by Theorem III.1, the incoherent state loses stability at Example III.3.(Classical Kuramoto model on full graph) In the case of the full graph (i.e., complete graph with uniform coupling strength), we have Clearly, A is a non-invertible operator, and the only eigenvalue of A is 1 (The eigenvalue equation A f = λ f implies that f must be a constant, say f 0 = 0, satisfying f 0 = λ f 0 .Thus, λ = 1.)Moreover, in the case that λ ∈ C \ {0, 1}, the operator A − λ I is invertible since for any g ∈ L 2 (Ω, m), the pre-image f is achieved under the unique choice Thus, we have σ (A) = {0, 1}.Hence, for the Kuramoto model on the full graph, we obtain by Example III.2 that in agreement with previous analysis (see also, e.g., Ref. 6).
Remark III.4.Sakaguchi 50 obtained for the critical coupling of the full graph the formula (in Sakaguchi's notation) In our framework, matching the assumptions and the notation correctly, we have ω 0 = 0, D = 1 β , and g = δ 0 .Note that in Sakaguchi's framework, the variance of the Brownian term f i (t) is 2Dt, while in our framework, the Brownian term 2β −1 W k has variance 2 β t for each k; thus, we must have D = 1 β .Thus, Sakaguchi's formula simplifies to which is exactly just the special case of the far more general formula we calculated in Example III.3.
Although we have now a very nice formula for C , it is not immediately clear for which classes of networks this formula works as N → ∞.After all, Theorem III.1 only makes claims about stability/instability based upon the assumption of the validity of the mean-field VFPE.Only if we already knew that the mean-field limit VFPE would be valid for certain classes of networks, i.e., if it does approximate -in a suitable sense -the oscillator system for finite but large N, then we could be certain applying our result for finite large networks.Proving such an approximation result in full generality is difficult, although first steps exist for the deterministic Vlasov case 43,44,46,47 .For example, one issue in this context is that the mean-field only holds in a scaling limit upon renormalizing the sums appearing in the Kuramoto model suitably via the density of the graph.However, empirically testing the formula for C via various classes of large finite networks using numerical simulation is certainly possible, and we shall proceed with this approach.

IV. INCOHERENCE-COHERENCE TRANSITION FOR FINITE AND INFINITE OSCILLATOR NETWORKS
We want to check the prediction for the incoherencecoherence transition given in Theorem III.1 for the mean-field limit by numerical simulations.The challenges we face in doing so stem from the fact that numerical simulations are bound to a finite-dimensional representation of Eqs. ( 1) and to a finite simulation time.Thus, while Theorem III.1 can only hold in an approximative sense for N < ∞, the finite system size and simulation time also incur uncertainty in the detection of the incoherence-coherence transition.Several points need to be taken into account when detecting the transition from incoherence to coherence that we outline below.
To see this, it is instructive to observe the dynamics for the case of the Kuramoto model where oscillators interact with D(u) = sin(u) on a complete graph with uniform coupling.The collective dynamics of all oscillators is described by the order parameter r(t) defined in (3) and is shown in Fig. 1 for numerical solutions of (2) for varying coupling strengths C and fixed system size N = 1000 and noise level β = 50.Initial phases are chosen to correspond to incoherent oscillations (see Sec. IV A on numerical methods).The dynamics of the order parameter r is subject to fluctuations, which stems from two sources: i) the stochastic dynamics inherent to the system and ii) finite-size effects induce pseudo-random fluctuations of order O(N −1/2 ) that vanish in the limit N → ∞ 51 .After a transient time, T tr , we observe that the dynamics settle into a quasi-stationary state (on average); i.e., the order parameter fluctuates around a constant mean value and is bounded by minimal and maximal values.If the trajectory after the transient attains a minimal value arbitrarily close to 0 during the observed time interval, we say that the population oscillates incoherently; if the minimal value never approaches 0, the dynamics are said to be (partially) coherent or synchronized (perfect synchrony occurs only for r = 1), and we observe increasing synchrony for larger C. Accordingly, Fig. 1 allows us to distinguish incoherent oscillations for weak coupling strengths (C = 0 to C = 0.03), and partially coherent oscillations occur for stronger coupling (C ≥ 0.04), which agrees well with the prediction of C = 0.04 given by ( 23) for the con-tinuum limit.For further details on the incoherence-coherence transition of the Kuramoto model, see also Ref. 6.
These observations point toward an implementation of numerical methods and measurements as outlined below.

A. Numerical methods
We calculate numerical solutions of (1) with a first-order Euler-Maruyama scheme with a time step ∆t = 0.01.Initial conditions/phases correspond to low synchrony compliant with incoherence, i.e., either the equidistant state θ k (0) := 2πk/N (Uniform complete graph, Erdös-Rényi graph, regular ring lattice with r = 400, spherical graph) or the random state where θ k (0) (Regular ring lattice with r = 25, sinuisodal graph, Lorentzian graph) is drawn from the uniform distribution on the interval [0, 2π) (two types of initial conditions were chosen since other attracting states were present for the regular ring lattice with r = 25).To characterize the post-transient dynamics, we use the order parameter r(t) = 1 N ∑ N j=1 e iθ j (t) in Eq. ( 3) and measure its temporal minimum and maximum, as well as its time average, where T := [T tr , T ] with T tr being the (estimated) transient time and T the total length of the simulation.
To average over stochastic effects, such as Brownian motion and random graphs (Erdös-Rényi and small-world), we average these measurements over several realizations of solutions of (2) (i.e., ten realizations to account for Brownian motion for eight (random) graph realizations) and denote ensemble averages with angular brackets • .To numerically test Theorem III.1, we calculate r min , r , and r max for different values of C and compare the resulting curves with C .The sampling points for the coupling C are non-uniformly spaced with a higher density in regions of interest (indicated as blue dots in Fig. 2).
A suitable transient time T tr can be determined based on the following considerations.The actual transient is maximal for C = C and decreases for C > C ; see Fig. 1.One could estimate T tr for each value of C individually to optimize for computational effort; but for simplicity, we estimated the length of T tr only at C = C and used this T tr for all probed values of C, as this choice guarantees a sufficiently long transient time.Due to the fluctuations present in the signal of r(t) (pseudo-random fluctuations and stochastic noise), the estimation of T tr is heuristic; i.e., it is done by visual inspection.This estimate of T tr improves with increasing N as the amplitude of (pseudo-random) fluctuations decreases.Taking these considerations into account, we chose T tr = 700, T = 1000, ∆t = 0.01, N = 1000 for all our numerical solutions of (2).

B. Graph(on) topologies and their associated incoherence-coherence transitions
We now define different graph structures for which we carry out numerical simulations to test for the onset of the incoherence-coherence transition.Results for the incoherence-coherence transitions for the various graph topologies are summarized in Fig. 2.

Incoherence-Coherence threshold for finite and infinite oscillator systems
We extend our analysis to different coupling topologies while using the coupling interaction D(u) = sin(u).The theoretically predicted threshold for the incoherence-coherence transition, C , valid in the mean-field limit is calculated using (22).We shall compare the numerical findings to this theoretical prediction for coupling topologies where it is possible, However, for certain graphops A, a characterization of σ (A) exceeds the scope of this study (Spherical graph in Sec.IV B 2 d; Lorentzian graph in Sec.IV B 2 f).In such cases, we instead compute the eigenvalues σ N (A N ) of a discrete coupling matrix A N that approximates A. We expect that the finite-dimensional matrices A N can be used to provide an approximation to (at least the boundary of) the spectrum of the limiting graphop A as N → ∞, and therefore, C can be approximated by its discrete corollary where Λ N (A N ) := max λ ∈σ N (N −1 A N ) |λ | is the maximal eigenvalue associated with A N .Finally, we also mention the possibility of "spectral pollution," 52 which, in principle, can occur when numerically approximating the spectrum of an operator with finite-dimensional matrices.However, as we shall see, our numerical and analytical results are consistent; we, therefore, anticipate that the numerical calculations are sufficiently stable.

Coupling topologies
a. Regular ring lattice with r neighbors.Nodes for this coupling topology may be imagined to be arranged on a ring such that every node is linked to a given number of r nearest neighbors.In the continuum limit N → ∞ , the ring lattice graphon can be defined as where 0 ≤ h ≤ 1/2 is the (continuous) coupling range for oscillators located at x and y on Ω.The graphop A defined via this graphon kernel K has Λ(A) = 2h (this can be shown, e.g., by writing K(x, y) as a Fourier series, and the values of σ (A) are given in Ref. 53.)For N < ∞, we simply define the regular ring lattice graph via where the (discrete) coupling range r = r(h) ∈ [N] for oscillators located at k and j in [N] satisfies 0 ≤ r ≤ N/2 with N even.It is easy to check that Λ N (A N ) = 2r/N.In our simulations, we choose r from which the value h for the corresponding graphon kernel follows via h = r/N.We then have Λ(A) = Λ N (A N ) = 2r/N.We note two limiting cases; namely, we obtain all-to-all coupling for h = 1/2 and zero coupling for h = 0. b.Erdös-Rényi graph.The Erdös-Rényi (ER) graph(on) is constructed in a random process where the presence (or absence) of every edge (of the complete graph) is chosen with a probability p ∈ [0, 1].
In the continuum limit, N → ∞, the Erdös-Rényi graphop simply becomes the complete (all-to-all) graphop with constant uniform coupling strength p; i.e., the corresponding graphon kernel is K(x, y) = p; see Ref. 54.It follows then that C = 2/(β p).
For finite oscillators N < ∞, a realization of the ER graph on N nodes may be obtained by drawing N(N − 1)/2 random numbers a k j , 1 ≤ k < j ≤ N, from the uniform distribution on the interval [0, 1].The adjacency matrix of the graph is then For p = 1, we obtain all-to-all coupling with uniform coupling strength 1 (complete graph), while p = 0 yields zero coupling.c.Small-world graph.The small world (SW) graph 55 interpolates between a regular ring lattice and a ER graph structure, thus creating a topology that is quite regular but also entertains random links across the network.This structure results in short path lengths even when nodes are far away on the ring.
For finite graphs, N → ∞, the small world graphop A can be constructed via the graphon kernel 56,57 given by where with (continuous) coupling range 0 ≤ h ≤ 1/2 (note that W (x, y) is identical to K(x, y) in ( 30) further above for the regular ring lattice).It can be shown that Λ(A) = 2h (to see this, one needs to write K(x, y) in terms of a Fourier series; the values of σ (A) are given by Gao and Caines 53 ).
For N < ∞, realizations of the SW graph on N nodes may be obtained via the procedure introduced by Watts and Strogatz 55 : One starts with a regular ring lattice on N nodes with r nearest neighbors (discrete coupling range).One selects a constant probability p ∈ [0, 1].For each node k and each link between k and its r nearest neighbors to the right, we draw a random number X ∈ [0, 1] i.i.d.from the uniform distribution.If X ≤ p, we draw a random integer j from the uniform distribution on [N].If k = j and the edge (k, j) does not yet exist, it is created and the old link deleted.
In our numerical setting, we simply pick a value r and the value h for the corresponding graphon kernel follows from h = r/N.We numerically confirmed that Λ N (A N ) ≈ Λ(A) = 2h.We shall thus use the value C ,N ≈ C = 1/(β h).
d. Spherical graph.The action of the spherical graphop A : L 2 (S 2 ) → L 2 (S 2 ) on a function f is defined by where ν x is the uniform measure.The spherical graphop thus integrates f over the circle on the unit sphere that consists of all the points perpendicular to x.The resulting circle is the equator of the point x.The spherical graphop does not have a graphon kernel or a known spectrum; therefore, we need to calculate C via (29).Moreover, a matrix approximation to the spherical graphop has to our knowledge not yet been proposed.Here, we propose a possible approximation without claiming any convergence properties as N → ∞.Choosing N (approximately equidistant) sample points x 1 , . . ., x N on the unit sphere, we may obtain a matrix A N approximating A by defining A N k j = A N jk = 1 if x k and x j are approximately perpendicular; otherwise, A N k j = A N jk = 0.The discretized version of (35) then reads We refer to A N as a spherical graph.Three requirements should be made on A N .For each point x k , the points x j for which A N k j = 1 should (i) lie sufficiently close to the equator of x k , (ii) be sufficiently equidistant, and (iii) be (almost) equally many for all k.Clearly, if we take an arbitrary point on the unit sphere, one can place M perfectly equidistant points on its equator.However, (i) and (ii) must be fulfilled reasonably well for all N points and their respective equators.Therefore, the points should form a regular grid.While a perfectly regular grid of N > 6 points on the sphere is impossible, there exist approximately regular grids 58 .Here, we place the points in a spiral of width 0.1 + 1.2N around the sphere, starting and ending (approximately) on the poles.This method is implemented in the Mathematica Software package 59 .We denote the set of points with this spacing on the unit sphere as P. The task is to determine subsets E k ⊂ P, 1 ≤ k ≤ N, such that each E k discretizes the equator of x k .To this end, we first calculate p k j := | x k , x j | for each 1 ≤ k < j ≤ N, to determine how close the pairs of points are to being perpendicular.Then, we specify M, the desired (approximate) cardinality of all E k 's.Now, we can, for each k, find the (approximately) M points x j with the smallest values of p k j and make them members of E k , under the constraint that if x k ∈ E j , then x j ∈ E k , to ensure that A N is symmetric.We end up with a (symmetric) A N that fulfills demands (i) and (ii) in an acceptable manner, while demand (iii) is fulfilled well: |E k | is either M or M − 1 for all 1 ≤ k ≤ N (see Fig. 3).
Since the spherical graphop cannot be defined via a graphon, we find Λ N (A N ) = 0.04998, and thus, C ,N = 0.8003.
e. Sinusoidal graph.In the sinusoidal coupling topology, nodes are coupled most strongly to their nearest neighbors, the coupling then smoothly decreases the farther neighbors are apart, finally the coupling is zero between nodes opposite on the ring.We define the graphon kernel as It can be shown 53 that the graphop A induced by K has Λ(A) = 1/2 so that C = 4/β = 0.08.For N < ∞, we define the matrix A N by f. "Lorentzian" graph.We also consider graphs for which a mean-field description is more challenging, and which, therefore, could potentially fail to exhibit the behavior predicted by Theorem III.1.A good candidate would be an irregular and sparse graph with few very strong links, while the vast majority of links is very weak.We can define such a topology based on the Lorentzian (graphon) kernel, where x 0 , y 0 ∈ [0, 1].K(x, y) peaks in the points (x 0 , y 0 ) and (y 0 , x 0 ), and converges to a sum of delta distributions centered at these points as µ → 0. We approximate this graphon in the finite representation as

V. CONCLUSION AND OUTLOOK
We formulated a mean-field theory for stochastic phase oscillator models with nontrivial coupling, i.e., heterogenous graph topologies and coupling weights.Our analysis for Kuramoto-type models with odd symmetric coupling functions, obtained via linearization around the incoherent solution branch, yields an exact formula for the critical coupling strength C at the incoherence-coherence transition in the mean-field limit.Numerically integrating finite representations (see Eq. ( 1)) agrees very well with the predicted threshold C (Eq. ( 22)) for a wide range of heterogeneous graph structures (see Fig. 2) 60 .We, therefore, expect our theory to be applicable to a large range of applications with heterogeneous oscillator interactions, such as systems with non-uniform coupling associated with chimera states 61 or XYoscillators type models with random coupling 34,62 .
For certain graph topologies characterized by strong sparsity, large variance in coupling strengths, or other types of "clusterization" implying coupling fragmentation in the network, the mean-field description is expected to break down, in particular, also in terms of correctly predicting the incoherence-coherence transition for finite-size systems.We found that such a problem occurs at least for one instance, namely, for the Lorentzian graph topology (see Fig. 2 (d) and 2(e)), for which the detection of a sharp transition point numerically is difficult.The Lorentzian graph is characterized by only a few nodes with very strong edge weights, while the vast majority of edge weights are very weak: the graph topology is effectively very sparse.This implies that we need a much larger C to observe coherent oscillations.As becomes apparent from comparing panels (d) and (e), the different quality of the incoherence-coherence transition between the Lorentzian and the other graphs considered is especially pronounced as the effective sparsity increases (µ → 0).Note that not merely larger overall coupling strength C is needed to achieve (partial) coherence, when compared to other topologies; if that were the case, one would just observe larger C for the Lorentzian graph as compared to the other graphs, and the coherence onset would still set in at C = C .Rather, the onset of coherence appears to be delayed beyond C = C so that the increase of partial synchrony sets in very slowly as C increases.This observation becomes especially pronounced for very small µ so that the coupling kernel becomes effectively very sparse.Thus, the Lorentzian graph represents an interesting coupling topology that demarcates a possible class of graphs for which -at least for certain values of µour mean-field description and prediction for the incoherencecoherence transition for the finite-size representation break down.
While we extended the mean-field theory for the stochastic Kuramoto model with all-to-all connectivity and uniform coupling strengths to heterogeneous connectivity with nonuniform coupling strengths, certain constraints apply to our model.These may limit the validity of our theory and prompt avenues for future research.For instance, we have assumed that the coupling function D(u) is odd.This assumption excludes, in particular, the Kuramoto-Sakaguchi model, which has a coupling function D(u) = sin (u + α) with a phase-lag α.This phase-lag allows one to tune the coupling interaction to be a sine function vs cosine, distinguishing gradient-like and integrable dynamics, respectively (compare with Eqs.(2) without noise (β −1 = 0) and implies different incoherence-coherence transitions (Note that a mix of such interaction is also essential to observe symmetry breaking chimera states with nonuniform synchronization patterns on the network 61,63 ) -extending our theory to such interactions would be of interest.Coupling functions D(u) of higher harmonic order have recently attracted much interest, which imply more complicated stability regimes and transitions between incoherence and coherence [64][65][66][67] .Moreover, interactions with arbitrary (e.g.non-symmetric) coupling interactions D = D(u k , u l ) are possible 68 which imply directed graph topologies 47 .While we studied the Kuramoto model with identical intrinsic frequencies, the presence of distributed frequencies is also of interest.Finally, extensions to other phase oscillator models, such as the Kuramoto model with inertia 26,69 or the theta neuron (or QIF neuron) that only performs rigid rotations corresponding to spiking above a threshold current, are worth mentioning.It would be very useful to derive rigorous mean-field descriptions for the above-mentioned systems; today, mean-field descriptions are available only for full graph structures 70,71 .Finally, one might also consider transitions between -or bifurcations of -states other than incoherence or coherence, such as chimera states or twisted states.Twisted states arise in bifurcations due to negative eigenvalues from the graph operator 72 .It would be interesting to extend the mean-field theory developed here to such cases.Some work in these directions has been done in the context hypergraphs 73 .
Another important avenue for future research is to clarify the validity regime for mean-field descriptions for very sparse and very heterogeneous structures.Note carefully that the effective dimension of the VFPE mean-field equation (15a) will grow the more heterogeneous the graph is due to the dependence of the node type encoded by points in Ω.Hence, a mean-field description can still exist, and our results indicate that this mean-field is often still very useful to determine whether some number of nodes starts to transition from incoherence to partial synchronization.Yet, for more complex patterns, involving an interplay between all different mean-field node types on very sparse structures, we anticipate that the mean-field description will eventually not be of much use as it is also high-dimensional.In summary, to fully determine the theoretical and practical limitations of heterogeneous meanfield VFPEs remains a challenging problem for future work.

VI. ACKNOWLEDGMENTS
MAG and CK gratefully thank the TUM International Graduate School of Science and Engineering (IGSSE) for support via the project "Synchronization in Co-Evolutionary Network Dynamics (SEND)."BJ and EAM acknowledge the DTU International Graduate School for support via the EU-COFUND project "Synchronization in Co-Evolutionary Net-work Dynamics (SEND)".CK also acknowledges partial support by a Lichtenberg Professorship funded by the Volkswagen Stiftung.

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

FIG. 1 .
FIG.1.Time evolution of r for numerical solutions of (2) for different values of C. Clearly, the time traces of r are subject to (random) fluctuations.Also, for higher C, the time traces of r settle, after a transient, at some quasi-stationary state (dashed lines).Other parameters are β = 50 and N = 1000.

FIG. 3 .
FIG. 3. Matrix approximation for the spherical graphop.Blue dots indicate the sample points, black circles mark points belonging to the discretized equator of one of the exemplary points (yellow diamond), and red squares mark points belonging to the discretized equator of the other exemplary point (green diamond).The two exemplary points are thus members of each other's discretized equators.Parameters are N = 1000 and M = 50.