The mean field analysis of the Kuramoto model on graphs II. Asymptotic stability of the incoherent state, center manifold reduction, and bifurcations

In our previous work [Chiba, Medvedev, arXiv:1612.06493], we initiated a mathematical investigation of the onset of synchronization in the Kuramoto model (KM) of coupled phase oscillators on convergent graph sequences. There, we derived and rigorously justified the mean field limit for the KM on graphs. Using linear stability analysis, we identified the critical values of the coupling strength, at which the incoherent state looses stability, thus, determining the onset of synchronization in this model. In the present paper, we study the corresponding bifurcation. Specifically, we show that similar to the original KM with all-to-all coupling, the onset of synchronization in the KM on graphs is realized via a pitchfork bifurcation. The formula for the stable branch of the bifurcating equilibria involves the principal eigenvalue and the corresponding eigenfunctions of the kernel operator defined by the limit of the graph sequence used in the model. This establishes an explicit link between the network structure and the onset of synchronization in the KM on graphs. The results of this work are illustrated with the bifurcation analysis of the KM on Erd\H{o}s-R{\' e}nyi, small-world, as well as certain weighted graphs on a circle.


Introduction
In 1970s, a prominent Japanese physicist Yoshiki Kuramoto described a remarkable effect in collective dynamics of large systems of coupled oscillators [13]. He studied all-to-all coupled phase oscillators with randomly distributed intrinsic frequencies, the model which now bears his name. When the strength of coupling is small, the phases are distributed approximately uniformly around a unit circle, forming an incoherent state. Kuramoto identified a critical value of the coupling strength, at which the the incoherent state looses stability giving rise to a stable partially synchronized state. To describe the bifurcation corresponding to the onset of synchronization, Kuramoto introduced the order parameter -a scalar function, which measures the degree of coherence in the system. He further showed that the order parameter undergoes a pitchfork bifurcation. Thus, the qualitative changes in the statistical behavior of a large system of coupled phase oscillators near the onset to synchronization can be described in the language of the bifurcation theory.
Kuramoto's discovery created a new area of research in nonlinear science [22,23,21]. The rigorous mathematical treatment of the pitchfork bifurcation in the KM was outlined in [5] and was presented with all details in [1]. The analysis in these papers is based on the generalized spectral theory for linear operators [2] and applies to the KM with intrinsic frequencies drawn from a distribution with analytic or rational probability density function. Under more general assumptions on the density, the onset of synchronization in the KM was analyzed in [7,9,6]. These papers use analytical methods for partial differential equations and build upon a recent breakthrough in the analysis of Landau damping [20].
In our previous work [3], we initiated a mathematical study of the onset of synchronization in the KM on graphs. Following [16,17], we considered the KM on convergent families of deterministic and random graphs, including Erdős-Rényi and small-world graphs among many other graphs that come up in applications. For this model, we derived and rigorously justified the mean field limit. The latter is a partial differential equation approximating the dynamics of the coupled oscillator model in the continuum limit as the number of oscillators grows to infinity. In [3], we performed a linear stability analysis of the incoherent state and identified the boundaries of stability. Importantly, we related the stability region of the incoherent state to the structural properties of the network through the spectral properties of the kernel operator defined by the limit of the underlying graph sequence [14]. In the present paper, we continue the study initiated in [3]. Here, we analyze the bifurcations at the critical values of the coupling strength K − c < 0 < K + c , where the incoherent state looses stability. As in the analysis of the original KM in [1], we have to deal with the fact that for K ∈ [K − c , K + c ] the linearized operator has continuous spectrum on the imaginary axis and no eigenvalues. Thus, neither the asymptotic stability of the incoherent state, nor the center manifold reduction for K = K ± c can be performed using standard methods of the bifurcation theory. To overcome this problem, following [1], we develop the generalized spectral theory, which allows effective analysis of the bifurcations at K ± c . This involves generalizing the resolvent operator and using it to define generalized eigenvalues. The generalized spectral theory is used i) to prove asymptotic stability of the incoherent state as a steady state solution of the linearized equation for K ∈ [K − c , K + c ]; ii) to obtain finite-dimensional center manifolds at the critical values of the coupling strength K = K ± c ; iii) to identify the bifurcations at K = K ± c . We do not develop nonlinear estimates in this paper. We also do not prove the existence of the center manifold rigorously. Both problems can be dealt with following the lines of the analysis in [1]. These very technical questions are left out to keep the length of this paper reasonable. Instead, we focus on the analysis of the bifurcations at K ± c . To this end, we develop all necessary tools for setting up the bifurcation problems, which are analyzed formally assuming the existence of the center manifold. The center manifold reduction yields stable spatial patterns arising from the bifurcations at the critical values of the coupling strength. We apply these results to the KM on several representative graphs, including small-world and Erdős-Rényi graphs. Our results for the KM on these and many other graphs are verified numerically in the follow-up paper [4]. In the remainder of this section, we review the model and the main outcomes of [3] and explain the main results of this work.
We begin with a brief explanation of the graph model that was used in [3] and will be used in this paper. In [3], we adapted the construction of W-random graphs [15] to define a convergent sequence of graphs. The flexible framework of W-random graphs allows us to deal with a broad class of networks that are of interest in applications. Let W be a symmetric measurable function on the unit square I 2 := [0, 1] 2 with values in [−1, 1] and let X n = {ξ n1 , ξ n2 , . . . , ξ nn } form a triangular array of points ξ ni , i ∈ [n] := {1, 2, . . . n}, n ∈ N, subject to the following condition (1.1) is a weighted graph on n nodes labeled by the integers from [n], whose edge set is Each edge {i, j} ∈ E(Γ n ) is supplied with the weight W nij := W (ξ ni , ξ nj ). In the theory of graph limits, W is called a graphon [14]. It defines the asymptotic properties of {Γ n } for large n.
Consider the initial value problem (IVP) for the KM on Γ ṅ The intrinsic frequencies ω i , i ∈ [n], are independent identically distributed random variables. The distribution of ω 1 has density g(ω). For the spectral analysis in Section 3, we need to impose the following assumptions on g: a) g : R → R + ∪ {0} is an even unimodal function, and b) g is real analytic function with finite moments of all orders: For instance, the density of the Gaussian distribution satisfies these conditions. The KM on weighted graphs {Γ n } (1.2), (1.3) can be used to approximate the KM on a variety of random graphs (cf. § 4.2 [3]).
Along with the discrete model (1.2) we consider the IVP for the following partial differential equation Here, ρ(t, θ, ω, x) is the conditional density of the random vector (θ, ω) given ω, and parametrized by (t, x) ∈ R + × I, and S = R/2πZ is a circle. In particular, interpreted as a probability measure on Borel sets A ∈ B(G), G = S × R × I, converges in the bounded Lipschitz distance [8] uniformly on bounded time intervals to the absolutely continuous measure provided µ n 0 and µ 0 are sufficiently close in the same distance. The latter can be achieved with the appropriate initial condition (1.3) and sufficiently large n (see [3,Corollary 2.3]). Therefore, the IVP (1.4),(1.5) approximates the IVP (1.2),(1.3) on finite time intervals for sufficiently large n.
An inspection of (1.4) shows that ρ u = 1/(2π), the density of the uniform distribution on S, is a steady state solution of (1.4). It corresponds to the incoherent (mixing) state of the KM. Numerics suggests that the incoherent state is stable for small K ≥ 0. The loss of stability of the incoherent state is interpreted as the onset of synchronization in the KM. This is the main focus of [3] and of the present paper. In [3], we identified the boundaries of the region of stability of the incoherent state in (1.4). Specifically, we showed , and is unstable otherwise.
The critical values K − c and K + c depend on the network topology through the eigenvalues of the compact symmetric operator W : L 2 (I) → L 2 (I) (1.10) The eigenvalues of W are real with the only accumulation point at 0. Denote the largest positive and smallest negative eigenvalues of W by µ max and µ min respectively. If all eigenvalues are nonnegative (nonpositive), we set µ min = −∞ (µ max = ∞). The main stability result of [3] yields explicit expressions for the transition points (1.11) Thus, the region of linear stability of ρ u depends explicitly on the spectral properties of the limiting graphon W . Recall that W represents the graph limit of {Γ n }. Thus, (1.11) links network topology to synchronization in (1.2). For the classical KM (all-to-all coupling), W (x, y) = 1 and µ min = −∞, µ max = 1, which recovers the known Kuramoto's formula.
In the present paper, we study the onset of synchronization in (1.2) in more detail. After some preliminaries and preparatory work in Sections 2 and 3, we revisit linear stability of the incoherent solution. This time, we show that despite the lack of eigenvalues with negative real part and the presence of the continuous spectrum on the imaginary axis, the incoherent state is an asymptotically stable solution of the linearized problem (cf. Theorem 4.1). This is a manifestation of the Landau damping in the KM.
In Section 5, we study the bifurcation at K + c with a one-dimensional center manifold. To this end, we recall the order parameter W (x, y)e iθ ρ(t, θ, ω, y)g(ω)dθdωdy, (1. 12) which was introduced in [3] as a measure of coherence in the KM on graphs. This is a continuous analog of the local order parameter 1 n n j=1 W nij e iθ nj (t) for the discrete model (1.2). The order parameter generalizes the original order parameter used by Kuramoto for the all-to-all coupled model. Note that (1.12) depends on x and contains information about the structure of the network through W . As will be clear below, the order parameter plays an important role in the analysis of the mean field equation. In particular, it can be used to locate nontrivial steady state solutions. To this end note that the velocity field (1.6) can be conveniently rewritten in terms of the order parameter (1. 13) In particular, for a given steady state of the order parameter written in the polar form (1.14) the velocity field takes the following form Setting ∂ θ (V ρ) = 0, we find the corresponding steady state solution of (1.4) where δ stands for the Dirac delta function. The stationary solution (1.15) has the following interpretation: the first line describes phase-locked oscillators, while the second line yields the distribution of the drifting oscillators. Thus, solutions of this form may combine phase-locked oscillators and those moving irregularly. Such solutions are called partially phase-locked or partially synchronized. The phase of a phase-locked oscillator at x with a natural frequency ω, is given by , (1. 16) provided |ω| ≤ KR(x). In Sections 5 and 6, we will identify branches of stable equilibria bifurcating from K ± c in terms of the corresponding values of the order parameter. Then Equation (1.16) will be used to describe the corresponding stable phase-locked solutions.
In Section 5, assuming that µ max is a simple eigenvalue of W, we show that the coupled system (1.2) undergoes a supercritical pitchfork bifurcation at K + c . Specifically, we derive an ordinary differential equation for the order parameter h and show that the trivial solution of this equation looses stability at K + c and gives rise to a stable branch of (nontrivial) equilibria, corresponding to partially synchronized state (cf. (5.1)). In Section 6, we consider the onset of synchronization in networks with certain symmetries (cf. (6.1)). This leads to the bifurcation with a two-dimensional center manifold. The bifurcation analysis in Sections 5 and 6 is illustrated with the analysis of the KM on Erdős-Rényi, small-world graphs, and to a class of weighted graphs on a circle.

Preliminaries
In the remainder of this paper, we will assume that K ≥ 0. The case of negative K is reduced to that above by switching to K := −K and W := −W . Furthermore, without loss of generality we assume that µ max > 0.

The eigenvalue problem
The multiplication operator M iω : H → H defined by is a closed operator. The continuous spectrum of M iω fills the imaginary axis Since P is compact (as a Hilbert-Schmidt operator), T : H → H is closed and σ c (T) = iR.
Next we turn to the eigenvalue problem where T and P are operators on H (cf. (2.10) and (2.5)).
We will locate the eigenvalues of T through the eigenvalues of W (cf. (1.10)). Since W is a compact symmetric operator on L 2 (I), it has a countable set of real eigenvalues with the only accumulation point at zero. All nonzero eigenvalues have finite multiplicity.
Suppose λ is an eigenvalue of T and v ∈ H is the corresponding eigenfunction. Then a simple calculation yields (cf. [3]) where Equation (2.14) yields the equation for eigenvalues of T where µ is a nonzero eigenvalue of W.
Using (2.17), we establish a one-to-one correspondence between the eigenvalues of W and those of T. Specifically, for every positive eigenvalue of W, µ, there is a branch of eigenvalues of T, Recall that µ max stands for the largest positive eigenvalue of W. Then for K ∈ [0, K(µ max )) there are no eigenvalues with positive real part. Furthermore, for small ε > 0 and K ∈ (K(µ max ), K(µ max ) + ε) there is a unique positive eigenvalue of T, λ(K, µ max ), which vanishes as K → K(µ max ) + 0 (see [3] for more details).

The generalized spectral theory
The major obstacle in studying stability and bifurcations of the incoherent state is the continuous spectrum of the linearized problem on the imaginary axis (cf. (2.12)). To deal with this difficulty, we develop the generalized spectral theory following the treatment of the classical KM in [1]. Below, we outline the key steps in the analysis of the generalized eigenvalue problem referring the interested reader to [1], [2] for missing proofs and further details.

The rigged Hilbert space
In this subsection, we define a rigged Hilbert space (a.k.a. Gelfand triple) [11] X ⊂ H ⊂ X , where H is a Hilbert space, and X is a dense subspace of H, whose topology is stronger than that of H. Throughout this paper, we assume that X is a locally convex Hausdorff topological vector space over C and X its dual space, the space of continuous antilinear functionals on X. Let ·, · denote the pairing between X and X, i.e., for l ∈ X and f ∈ X, l, f := l(f ) stands for the corresponding antilinear functional. To use the generalized spectral theory (cf. [2]) we also need X to be a quasi-complete barreled space.
We take L 2 (R × I, g(ω)dωdx) as the Hilbert space H and contruct X as follows. Let Exp(β, n) be the set of holomorphic functions on the region C n := {z ∈ C | Im(z) ≥ −1/n} such that the norm ||φ|| β,n := sup is an increasing sequence in n. By Montel's theorem, the inclusion Exp(β, n) → Exp(β, n + 1) is a compact operator for any n ≥ 1. By Komatsu's theorem [12], the inductive limit is a complete Montel space. In particular, it is a complete barreled (DF) space. Similarly, the inductive limit Exp(β) is a complete barreled (DF) space. The properties of Exp are described in detail in [1].
Let Exp(β, n) ⊗ L 2 (I) be a projective tensor product. Since the identity map L 2 (I) → L 2 (I) is weakly compact, the inclusions Exp(β, n)⊗L 2 (I) → Exp(β, n+1)⊗L 2 (I) and Exp(β)⊗L 2 (I) → Exp(β +1)⊗ L 2 (I) are weakly compact operators. By Komatsu's theorem [12], the inductive limit X := Exp ⊗ L 2 (I) is a complete barreled (DF) space and X is a Fréchet space. For every f ∈ X we have f (ω, ·) ∈ L 2 (I) for each ω ∈ R. In addition, f (ω, x) is holomorphic in ω on the upper half plane, where it can grow at most exponentially. Then the operator T and the rigged Hilbert space X ⊂ H ⊂ X satisfy all assumptions of the generalized spectral theory in [2].
. Thus, l ∈ H can be viewed as an element of X .

The generalized eigenvalue problem
In this subsection we calculate the resolvent of T and spectral projections. With the rigged Hilbert space defined above, we will view the resolvent as an operator from X to X .
Below, we will need to construct analytic continuation for certain functions involving integrals of Cauchy type. For this, we are going to use an implication of the Sokhotski-Plemelj formulas, which we formulate as a separate statement for convenience.
is an analytic function in the right and left open half-planes of C. Furthermore, for z = x+iy, the following formulas determine the limits of F (z) as x → 0±:

4)
is an entire function.

The generalized resolvent
Our next goal is to compute the resolvent of T R(λ) = (λ − T) −1 . (3.5) To this end, we first compute R(λ) for Re(λ) > 0 and extend it analytically to the left half-plane as an operator from X to X .
In the right half-plane Re(λ) > 0, R(λ) can be rewritten as follows where and I stands for the identity operator. Note that A(λ) ceases to exist as the multiplication operator on H as Re(λ) → 0 (recall that the imaginary axis is the continuous spectrum of M iω ). However, it can be extended to the left half-plane as as an operator A : X → X defined as follows where P × : X → X is the dual operator of P.
Since T has the continuous spectrum on the imaginary axis, R(λ) can not be continued to the left-half plane as an operator on H. We define the generalized eigenvalues of T as the singularities of the generalized resolvent R(λ).

10)
In this case, v is called a generalized eigenfunction.
Remark 3.4. Since the range of the operator P × : X → X is in X, A(λ)P × v is well-defined for v ∈ X .   where µ is a nonzero eigenvalue of W and Re(λ) = 0, (3.12) The right hand side of (3.12) is an entire function (cf. Corollary 3.2). For Re(λ) > 0, (3.11) is reduced to the equation for the eigenvalues of T (cf. (2.17)). In this case, the corresponding generalized eigenfunction v is included in L 2 (R × I, gdωdx), i.e., λ is an eigenvalue of T. On the other hand, for Re(λ) ≤ 0, the generalized eigenfunction v is not in H but is an element of the dual space X .
Since the generalized eigenvalue of T, λ, is a root of (3.11), R(λ)u is an X -valued meromorphic function for each u ∈ X. For Re λ > 0, it coincides with the restriction of R(λ) to X. Thus, R(λ) is a meromorphic continuation of R(λ) from the right half-plane to the left-half plane as an X -valued operator.

The generalized Riesz projection
Let µ be a positive eigenvalue of W and w ∈ L 2 (I) be the corresponding eigenfunction. The largest positive eigenvalue of W and the corresponding eigenfunction are denoted by µ max and w max respectively. For every K > K + c = 2/(πg(0)µ max ) there is a real positive eigenvalue of T, λ = λ(µ, K). The corresponding eigenfunction is given by As K approaches the critical value K + c from above, the eigenvalue λ(µ max , K) converges to 0+ along the real axis and at K = K + c it hits the continuous spectrum on the imaginary axis. The corresponding eigenfunction approaches the critical vector where lim* stands for the limit in X with respect to the weak dual topology 1 , i.e., the action of v + c ∈ X on u ∈ X is given by (3.15) 1 {ln} ⊂ X converges to l ∈ X if ln, f ∈ C tends to l, f for every f ∈ X.
Let λ ∈ C be a generalized eigenvalue of T. Then the generalized Riesz projection Π λ : X → X is defined by where γ(λ) is a simple closed curve around λ oriented counterclockwise that does not encircle or intersect the rest of the spectrum. Below, we shall refer to such curves as contours. The image of Π λ gives the generalized eigenspace of λ [1].
Theorem 3.8. Suppose the algebraic and geometric multiplicities of µ max coincide. Then the generalized Riesz projection of λ = 0, the generalized eigenvalue of T for K = K + c , has the following form

17)
where g 1 = − lim λ→0+ D (λ) −1 is a positive constant, A(λ) was defined in (3.7), andΠ µ stands for the Riesz projection onto the eigenspace of W corresponding to the eigenvalue µ. The operator D(λ) on H is defined by The proof of Theorem 3.8 relies on three technical lemmas. Below we state and prove these lemmas first and then prove the theorem.
Lemma 3.10. Let λ = λ(µ, K) > 0 be an eigenvalue of T corresponding to the positive eigenvalue of W, µ, and K > K + c , and suppose that the geometric and algebraic multiplicities of µ coincide. Then
Since the algebraic and geometric multiplicities of µ are equal, the singularity of (ζ − W) −1 at ζ = µ is a simple pole, and the other factor in the integrand of the above is regular at ζ = µ. Therefore, the right-hand side of (3.27) simplifies to By multiplying both sides of (3.28) by (2πi) −1 , we have Finally, sinceΠ µ is the projection on the eigensubspace of W, Thus, Proof. Using integration by parts n times, we obtain The application of Lemma 3.1 to the integral on the right-hand side yields (3.30).
Below will need the following implications of Lemma 3.11.
Proof. Differentiating D(z) and using (3.30), for z off the imaginary axis we have The integral on the right-hand side is of Cauchy type and Lemma 3.1 applies. By (3.3), Since g is even, g (0) = 0 and g is odd. Because g is also nonnegative and unimodal g (x) ≤ 0, x > 0. Thus,

Asymptotic stability of the incoherent state
We now return to the problem of stability of the incoherent state. Recall that in the Fourier space the incoherent state corresponds to the trivial solution Z = (z 1 , z 2 , · · · ) = 0. The linearization about Z = 0 shows that it is a neutrally stable equilibrium of (2.8), (2.9) for 0 ≤ K < K + c . There are no eigenvalues of T for these values of K and the continuous spectrum fills out the imaginary axis. Nonetheless, we show that the incoherent state is asymptotically stable with respect to the weak dual topology.
Theorem 4.1. For K ∈ [0, K + c ) the trivial solution of (2.8), (2.9) is an asymptotically stable equilibrium for initial data from X ⊂ X with respect to the weak dual topology on X .
Remark 4.2. The stability with respect to the weak dual topology is weaker than that with respect to the topology of the Hilbert space H. Still it is a natural topology for the problem at hand. In particular, Theorem 4.1 implies that the order parameter evaluated on the trajectories of the linearized problem tends to 0 as t → ∞.
By the Riemann-Lebesgue lemma, (4.1) We now turn to (2.8). By the Hille-Yosida theory, operator T generates a C 0 -semigroup e Tt , which can be computed using inverse Laplace transform (cf. where a > 0 is arbitrary. Thus, the (continuous) spectrum of T lies to the left of the integration path along x = a (see Fig. 1a).
For arbitrary φ, ψ ∈ H, we have For φ, ψ ∈ X, (λ − T) −1 φ, ψ H is an analytic function in the right half-plane, which can be extended to the entire complex plane as a meromorphic function R(λ)φ, ψ . Thus, Let K ∈ [0, K + c ) be fixed. Next we claim that one can choose ε = ε(K) > 0 such that there are no generalized eigenvalues of T on or inside the contour (Fig. 1b) for every R > 0. To construct C ε,R with the desired property, we first fix δ > 0. Then we recall that generalized eigenvalues of T satisfy (3.11). From (3.12), under our assumptions on g, there exists R 0 = R 0 (δ) such that there are no roots of (3.11) in the region because (3.11) can be reduced to 2/(Kµ) = O(1/|λ|) in D + R,δ for R 1. On the other hand, D(λ) is holomorphic. Thus, the set of roots of (3.11) (i.e., the set of generalized eigenvalues) does not have accumulation points in D − R 0 ,δ = {z ∈ C : |z| ≤ R 0 & − δ < Re(z) ≤ a}. Thus, we can choose ε > 0 such that there are no generalized eigenvalues in D + R 0 ,ε ∪ D − R 0 ,ε . This completes the construction of C ε,R with the desired property for any R > 0.
By the Cauchy Integral theorem, for any R > 0, and a+iR a−iR The integral on the left-hand side of (4.6) exists, by the Hille-Yosida theory. Therefore, the integrals on the right-hand exist too. Below, we show that the last two integrals on the right-hand side of (4.6) tend to 0 as R → ∞. Sending R → ∞ in (4.6) and using (4.4), we arrive at as t → ∞ because the integral exists and also tends to zero as t → ∞ due to the Riemann-Lebesgue lemma.
It remains to prove the following lemma.
Proof. We show that the integral −ε+iR a+iR e λt R(λ)φ, ψ dλ tends to zero as R → ∞. The second integral a−iR −ε−iR e λt R(λ)φ, ψ dλ can be treated in the same way. Further, we decompose the integral into two integrals as We show that the first integral on the right hand side tends to zero as R → ∞. For Re(λ) > 0, we have see (3.19). For the first term, we have iR a+iR e λt A(λ)φ, ψ dλ Since the integral above is finite, for any ε 0 > 0, there exists L > 0 such that 0 a e λt I |ω|>L On the other hand, the integrand as R → ∞ uniformly in x ∈ I, ω ∈ (−L, L) and λ ∈ (0, a). This implies that the integral iR a+iR e λt A(λ)φ, ψ dλ → 0, as R → ∞.
The singularity of φ λ is a generalized eigenvalue of T (cf. (2.14)). For 0 < K < K + c , there are no generalized eigenvalues of T in the right half-plane and on the imaginary axis. Further, D(λ) → 0 as |λ| → ∞ and, thus, D(λ) → 0 too. This shows that φ λ is bounded uniformly in λ on the region Re(λ) ≥ 0. By replacing φ with φ λ in the first estimate of the integral of A(λ)φ, ψ , we find that iR a+iR e λt A(λ) φ λ , ψ dλ tends to zero. This shows that iR a+iR e λt R(λ)φ, ψ dλ decays to zero as R → ∞. The second integral in (4.9) is analyzed in similarly. This completes the proof of Lemma 4.3.

Bifurcation with a one-dimensional null space
In the previous section, we proved asymptotic stability of the equilibrium at the origin of the linearized system (2.8), (2.9) for K ∈ [0, K + c ). On the other hand, for K > K + c there is a positive eigenvalue in spectrum of the linearized problem (cf. [3]). This signals a bifurcation at K + c . This bifurcation is analyzed in this present section. As in the classical KM, the loss of stability of the incoherent state at K + c and the development of partial synchronization for K > K + c is best seen in terms of the order parameter.
Throughout this section, we assume that the largest positive eigenvalue µ max of W with the eigenfunction w max is simple. Furthermore, we assume that at K + c there is a (one-dimensional) smooth center manifold of the equilibrium at the origin of (2.6), (2.7) 2 . Under these assumptions, below we show that the order parameter undergoes a supercritical pitchfork bifurcation at K + c . The stable branch of equilibria bifurcating from 0 is given by Formula (5.1) generalizes the classical Kuramoto's formula describing the pitchfork bifurcation in the allto-all coupled model to the KM on graphs. The network structure enters into the description of the pitchfork bifurcation through the largest eigenvalue µ max and the corresponding eigenspace.

Preparation
Throughout this section, we assume that µ max is a simple eigenvalue of W. Let K = K + c + with 0 < 1 and rewrite (2.6),(2.7) as follows where T 0 is T evaluated at K = K c and T = T 0 + P/2.
For small > 0, the equilibrium of (5.3), (5.4) at the origin has a 1D unstable manifold. We reduce the dynamics on the 1D unstable manifold, which we approximate by the center manifold of the origin for K = K + c , i.e., for = 0. For the latter, we assume z k = h k (z 1 ), k = 2, 3, . . . , on the center manifold, where h k are smooth functions such that h k (0) = h k (0) = 0.
Let Π 0 be the projection to the eigenspace of λ = 0 spanned by v + c (cf. Section 3.4). To track the evolution on the slow manifold we adopt the following Ansatz: where α > 0 is a small parameter, c(t) is the coordinate along the center manifold, and v + c is the generalized eigenfunction of T 0 corresponding to the zero eigenvalue (cf. (3.14)). The Ansatz (5.5)-(5.7) follows right away once existence of the center manifold is shown.

The slow manifold reduction
Projecting both sides of (5.3) onto the center subspace, we have To evaluate (5.17), we take the following steps Finally, Similarly, to evaluate we first compute It is instructive to recast (5.23) in terms of the order parameter h (cf. (2.3)). By Lemma 5.1, Thus, by multiplying both sides of (5.23) by αw max and neglecting higher order terms, we obtaiṅ (5.24) Equation (5.24) shows that the trivial solution (the incoherent state) looses stability at = 0 and for small > 0 there is a nonzero stable equilibrium (5.25)

Examples
In [3], we derived the transition formulas for the onset of synchronization in the KM on several networks. We now return to these examples and describe the transition to synchronization in more detail using the results of this section.
We start with the KM on the Erdős-Rényi graphs. To this end, let W ≡ p ∈ (0, 1). In [3], we showed that the largest positive eigenvalue of W in this case is µ max = p. The corresponding eigenfunction w max is constant. This yields the critical value K + c = 2(πg(0)p) −1 . By plugging in these values into (5.1), we obtain We next turn to the KM on small-world graphs. This family of graphs is defined via the following graphon: where p, r ∈ (0, 1/2) are two parameters. The former stands for the probability of long range random connections and the latter is the range of regular local connections (cf. [18]).
The largest eigenvalue of W p,r is equal to 2r + 2p − 4pr and the corresponding eigenfunction w max is constant (cf. [3]). This implies that the critical value is .
Using (5.1), we further have Bifurcation with a two-dimensional null space 6.

The slow manifold reduction
Many networks in applications can be described with the limiting graphon of the following form . The graphons of this form are used in the description of the small-world and many other networks (cf. §5.3 [3]).
A graphon satisfying (6.1) admits Fourier series expansion By Parseval's identity, It follows from (6.2) that the eigenvalues of the kernel operator W coincide with the Fourier coefficients c k , k ∈ Z. The Fourier modes e 2πikx , k ∈ Z, yield the corresponding eigenfunctions.
We continue to assume that the largest eigenvalue of W is positive, i.e., there is at least one positive coefficient c k . In view of (6.3), there is a finite set If M = {0} the null space of W is one-dimensional. This case was analyzed in the previous section. Here, we assume |M | = 2, i.e., there exists a unique m ∈ N such that µ max = sup{c k : k ∈ Z} = c m = c −m .
The corresponding eigenspace is spanned by w + = e 2πimx and w − = e −2πimx .
The linearization of (6.18) about these fixed points yields respectively. Thus, the second and the third fixed points are stable for 0 < K − K + c 1. This proves that the order parameter (6.9) tends to as t → ∞, where φ is a constant which depends on an initial condition.
The analysis of this section then yields two stable branches of solutions bifurcating from h ≡ 0 at K + c = 4(πg(0)) −1 ≈ 3. where the phase shift φ is determined from the initial condition. For small κ > 0, the system has a family of stable partially phase-locked solutions (1.15), which can be described as follows. The oscillators split into two groups depending on their intrinsic frequencies. If |ω i | < K κ p 2 , the oscillator i approaches one of the two phase-locked solutions: where Y (ω i ) = arcsin ω i K p 2 κ ∈ (−π/2, π/2) is a function of the random intrinsic frequency ω i . The oscillators in this group form a noisy 1-twisted state [19]. The oscillators with intrinsic frequencies |ω| > K κ p 2 are randomly distributed around S. The density of this distribution is given in the second line of (1.15). Figure 2 presents results of numerical integration of the KM with graphon (6.21) and randomly distributed intrinsic frequencies. The plots in Figure 2 a-c show asymptotic states of the KM for three increasing values of K, starting with K = 3.5 just near the critical value K + c ≈ 3.2. In Figure 2a there are many oscillators spread around S. However, the group of oscillators concentrating about a 1-twisted state is already visible. For larger values of κ, the twisted state becomes more pronounced (see Figure 2b,c). Twisted states bifurcating from the incoherent state are also present in the KM on small-world graphs (see [4] for the analysis of the small-world network and other examples).