Stability and bifurcation for the Kuramoto model

We study the mean-field limit of the Kuramoto model of globally coupled oscillators. By studying the evolution in Fourier space and understanding the domain of dependence, we show a global stability result. Moreover, we can identify function norms to show damping of the order parameter for velocity distributions and perturbations in $\mathcal{W}^{n,1}$ for $n>1$. Finally, for sufficiently regular velocity distributions we can identify exponential decay in the stable case and otherwise identify finitely many eigenmodes. For these eigenmodes we can show a center-unstable manifold reduction, which gives a rigorous tool to obtain the bifurcation behaviour. The damping is similar to Landau damping for the Vlasov equation.

The finite model consists of N oscillators with natural frequencies (ω i ) N i=1 whose phase (θ i ) N i=1 evolves as for i = 1, . . . , N , where K is the coupling constant and is the order parameter and the main observable of the system. This has the structure of a kinetic particle system where the phase angle θ corresponds to the spatial position and the natural frequency ω to the velocity. Following this structure, we will use the terminology of spatial and velocity distribution.
We study the mean-field limit for N → ∞, cf. [24,35]. Main features of the evolution were already identified in early studies [22,23] and [34] contains a more precise review. The spatially homogeneous distribution is a trivial stationary solution with η = 0 and is also called the (fully) incoherent state. The system then has a critical coupling constant K c so that for K < K c the incoherent state appears to be stable and for larger K the system bifurcates to partially locked states.
Even though the order parameter suggests a stable behaviour for K < K c , the linearised evolution operator in normal physical function spaces has a continuous spectrum along the imaginary axis and thus does not show decay directly [35]. Already in 1992 it was realised that the behaviour is similar to Landau damping in plasma physics [36].
Here our starting point is to study the model in Fourier variables for the phase θ and the frequency ω. By understanding the domain of dependence, we can prove a global stability result: There exists K ec such that for couplings K < K ec the evolution is controlled by ∞ 0 |η(t)| 2 dt < ∞.
For general velocity distributions K ec < K c , but for the Lorentzian and Gaussian distribution K ec = K c . As far as the author is aware, this is the first global result controlling the order parameter.
For the local behaviour of the spatially homogeneous state, we first note that the theory of the Volterra equation [19] gives a rigorous stability result for the linearised evolution of the order parameter η. In particular, it shows exponential decay when the Fourier transform decays exponentially and algebraic decay t −n when the Fourier transform decays like t −n . In particular the condition for algebraic decay is satisfied if the velocity distribution and the perturbation are in the Sobolev space W n,1 .
By working in Fourier variables, we can now define norms for the perturbation in which the linearised evolution has the same decays as the order parameter. Moreover, for exponential decay or algebraic decay t −n for n > 1, we are able to extend the linear stability to nonlinear stability of perturbations. In particular, this shows that the order parameter decays as predicted by the linear analysis.
In case the linear stability criterion fails and we have exponential regularity, the linear analysis allows us to identify finitely many eigenmodes. By constructing suitable norms we are now able to prove a center-unstable manifold reduction of the dynamics. In particular, this allows the rigorous derivation of the bifurcation and shows stability of the partially locked states in the supercritical bifurcation.
We stress that our results and used norms are independent of the particular velocity distribution. Thus the stability result and the center-unstable manifold reduction hold for all sufficiently regular distributions. To the best of the author's knowledge this result was unknown.
For Gaussian and rational velocity distributions a similar stability and bifurcation result has been obtained by Chiba [10,12]. The dynamics of the system can also be reduced by using a symmetry reduction, the so called Ott-Antonsen ansatz [28][29][30]. For special rational velocity distribution ((1+ω 2 ) −1 , (1+ω 4 ) −1 ) [28] and for the Bi-Cauchy distribution [25] this has been used to obtain the dynamics. In case of identical oscillators (i.e. the frequency distribution is concentrated at one frequency), it has been proved that the system will asymptotically be fully synchronised [20]. This result has been extended to compactly supported velocity distributions with sufficiently coherent initial data and large enough coupling constant [9].
The nonlinear stability and bifurcation have been studied in the modified system with white noise [32]. In this case the mean-field limit becomes a Fokker-Planck equation and the results have been obtained in the 1990s [7,15,35]. The concrete bifurcation behaviour for different velocity distributions and different interactions has been studied in more detail [1,6,13,14], see [2] for a review.
The above results show that the Kuramoto model is an interesting model to understand Landau damping where we can obtain more results than for the classical Vlasov equation for plasmas [4,26]. Other models are the 2d Euler equation [3] and the Vlasov-HMF model, where recently stability has been proved with Sobolev norms [17].
Remark. After finishing the work, the author got aware of the independent work by Fernandez, Gérard-Varet, and Giacomin [18]. This very interesting work shows the nonlinear stability for perturbations in the class C n (n ≥ 4) and follows more closely [17]. After finishing a first draft another preprint [5] appeared showing for analytic initial data exponential convergence to the incoherent state if the coupling constant K is small enough.
1.2. Strategy and detailed results. For finitely many oscillators with phase angles (θ i ) N i=1 and frequencies (ω i ) N i=1 we can look at the empirical measure ρ N,t (θ, ω) = N −1 N i=1 δ θ=θi(t) δ ω=ωi(t) in the phase space Γ = T×R, where T is the torus [0, 2π) for the phase angles. The set of probability measures on a space X is denoted by M(X) and in the mean-field limit N → ∞ we suppose that the initial distribution converges weakly to a probability distribution ρ in ∈ M(Γ). Then for all later times t ∈ R + the distribution converges weakly to a probability distribution ρ t and satisfies a PDE in a weak sense. For the limiting case N → ∞ we can thus study the family ρ t whose evolution is uniquely determined by the PDE. Given initial data ρ in , we call this family ρ t the solution to the Kuramoto equation or the mean-field solution.
Assuming a density function f (t, ·, ·) of ρ t , the system evolves formally as In the general case, we rigorously understand this PDE in a weak sense for measures. As a Cauchy problem, we will always assume that the initial data ρ in are given at time t = 0 and then the solutions ρ · will be defined for all t ∈ R + .
The velocity distribution g is a probability distribution M(R) which is constant in time, so that g(A) = ρ t (T × A) holds for all A ∈ B(R) and time t, where B(R) denotes the Borel sets of R.
Working formally with a density f , the velocity distribution has a density g and we can take the Fourier series for the spatial modeŝ for l ∈ Z. These then evolve as Since f is non-negative and its marginal g(ω) = θ∈T f (t, θ, ω)dθ is constant in time, we find for all time t f (t, 0, ω) = g(ω) and |f (t, l, ω)| ≤ g(ω) for all l ∈ Z.
Moreover, since f is real,f (l, ω) =f (−l, ω) so that it suffices to consider l ∈ N. The spatially homogeneous distribution is characterised by vanishing non-zero moments, so that the stability of the spatially homogeneous distribution is equivalent to the decay off restricted to l ≥ 1. The evolution of the restriction satisfies Exploiting that the non-linear interaction is skew-Hermitian, we find With K > 0 the system therefore is not dissipative and the convergence to zero must happen in a weaker sense. Moreover, it shows that the model is not timereversible, as K needs to change sign as well. Finally, for negative K this already shows global stability as then |η(t)| 2 dt is bounded. Hence the interesting case is K ≥ 0 which we assume in the following. In order to understand the evolution of the spatial modes l ≥ 1, we also Fourier transform the velocity variable ω as The difference to the spatially homogeneous state is again captured by the restriction to l ≥ 1. Hence the evolution of u(t, ·, ·) ∈ C(N × R) describes the evolution of the perturbation and is determined by (3) In this formulation the free transport always moves ξ → u(t, l, ξ) to the left and the interaction at (t, l, ξ) is given by u(t, 1, 0) = η(t) and u(t, l ± 1, ξ). Hence with M ≥ 0, the region N×{ξ : ξ ≥ −M } is its own domain of dependence. This already suggests to use norms focusing on ξ ≥ 0 in order to show convergence.
For the global energy estimate, we note that the interaction is skew-Hermitian except for the first mode. Thus we consider the energy functional where φ is an increasing weight. By balancing the gain of the first linear interaction term with the decay by the increasing weight, we can show our first main result.
Theorem 1. For a velocity distribution g ∈ M(R) suppose ∞ 0 |ĝ(ξ)|dξ < ∞ and let If the coupling constant K satisfies K < K ec , then there exists a finite constant c > 0 and a bounded increasing weight φ ∈ C 1 (R + ) with φ(0) = 1 so that for a solution to the Kuramoto equation with velocity marginal g the energy functional I(t) defined by Equation (4) satisfies for all t ∈ R + . In particular this shows that I is non-increasing and the order parameter η satisfies If, moreover, the initial distribution ρ in has a density f in ∈ L 2 (Γ), then . In order to quantify the decay and stability we introduce, for a non-negative weight φ, the weighted spaces L 1 (X, φ) and L ∞ (X, φ) for functions f : X → C with finite norm respectively. As weights we will use exp a defined as exp a (x) = e ax and p A,b defined as p A,b (x) = (A + x) b with the special case p b = p 1,b . In Fourier variables the linear evolution is determined by the PDE Here all modes decouple and only the first mode differs through the interaction. We denote its evolution operator by e tL , i.e. t, l, ξ → (e tL u in )(l, ξ) is the weak solution to Equation (5) in C(R + × N × R) with initial data u in given at time t = 0. As in the case of the Vlasov equation [31] the order parameter under the linear evolution is given by a Volterra equation.
Lemma 2. Letĝ ∈ C(R). Then the evolution operator e tL is well-defined from continuous initial data to the unique continuous solution. For u in ∈ C(N × R) let ν(t) = (e tL u in )(1, 0). Then ν satisfies the Volterra equation with the convolution kernel and and for l ≥ 2 (e tL u in )(l, ξ) = u in (l, ξ + tl).
Here the convolution f * g between two functions f, g ∈ L 1 The solution to the Volterra equation is given by the resolvent kernel r ∈ L 1 loc (R) as Hence the order parameter ν under the linear evolution can decay if u in (1, ξ) is decaying as ξ increases. We can express this for a ≥ 0 as and for A ≥ 1 and b ≥ 0 as With the same decay inĝ we have a sharp criterion for the resolvent kernel by the Paley-Wiener theorem for Volterra equations [ Here Lk denotes the Laplace transform of k, i.e.
Hence the linear evolution of the order parameter is stable if (Lk)(z) = −1 for all z ∈ C with z ≥ 0. By bounding the absolute value, we see that if K < K ec the system is linearly stable. In case of the Gaussian or the Cauchy distribution,ĝ is always positive so that K ec equals the critical coupling. Like in the Vlasov equation [31], this condition can be related to a Penrose criterion for the complex boundary, which also visualises how often the condition fails. Under this we can also relate the stability to the known condition [35,36].
For small perturbations, our second main result is the nonlinear stability of the incoherent state. In a first version we propagate control in sup l≥1 sup ξ∈R |u(t, l, ξ)|e a(ξ+tl/2) . Theorem 4. Let g ∈ M(R) and a > 0 such thatĝ satisfies the exponential stability condition of rate a in Theorem 3, i.e.ĝ ∈ L 1 (R + , exp a ) and (Lk)(z) = −1 for all z ∈ C with z ≥ −a. Then there exists a δ > 0 and finite c 1 such that for initial data ρ in ∈ M(Γ) with velocity marginal g, Fourier transform u in and M in := sup l∈N sup ξ∈R |u in (l, ξ)|e aξ ≤ δ the Fourier transform u of the solution to the Kuramoto equation with initial data ρ in satisfies sup t∈R + sup l≥1 sup ξ∈R |u(t, l, ξ)|e a(ξ+tl/2) ≤ (1 + c 1 )M in .
In particular, this shows that the order parameter η(t) = u(t, 1, 0) decays as O(e −at ).
For the algebraic decay we propagate control in Theorem 5. Let g ∈ M(R) and b > 1 such thatĝ satisfies the algebraic stability condition of parameter b in Theorem 3, i.e.ĝ ∈ L 1 (R + , p b ) and (Lk)(z) = −1 for all z ∈ C with z ≥ 0. For α > 0 satisfying b > 1 + α, there exists δ > 0 and finite γ 1 such that for initial data ρ in ∈ M(Γ) with velocity marginal g, Fourier transform u in and M in := sup l∈N sup ξ∈R |u in (l, ξ)|(1 + ξ) b ≤ δ the Fourier transform u of the solution to the Kuramoto equation with initial data ρ in satisfies In particular, this shows that the order parameter η(t) = u(t, 1, 0) decays as O(t −b ).
Note that the exponential and algebraic control are one-sided and thus do not control the L 2 or Sobolev norm in the original space. This explains why we can have decay with these norms, even though there is no decay of perturbations in the L 2 norm and the linear system is only neutrally stable [35]. For the application, however, the Sobolev norm W b,1 controls M in and g ∈ W b,1 ensures the required decay ofĝ.
If the stability condition fails, we can identify growing modes in the Volterra equation.
Theorem 6. For a ∈ R supposeĝ ∈ L 1 (R + , exp a ) and that the convolution kernel k defined by Equation (7) satisfies (Lk)(z) = −1 for z ∈ C with z = −a. Then the solution to the Volterra equation (6) is given by where r s is vanishing for negative arguments and satisfies r s L 1 (R + ,exp a ) < ∞ and λ 1 , . . . , λ n are the finitely many zeros of 1 + (Lk)(z) in z > −a with multiplicities p 1 , . . . , p n and b i,j depends on the residues of [1 + (Lk)(z)] −1 at λ i . This says that the linear evolution is governed by eigenmodes t j e λit for i = 1, . . . , n and j = 0, . . . , p i − 1 and a remaining stable part r s . The previous linear stability theorem is included in this formulation as r = r s .
Ifĝ ∈ L 1 (R + , exp a ) for a > 0, this shows that if the stability condition fails, there exists an unstable mode. Without extra assumption, for every a < 0 we havê g ∈ L 1 (R + , exp a ) because g is a probability distribution. Hence, unless we are at the critical case, we have unstable modes if the stability condition fails. Therefore, the linear stability condition is sharp.
The additional restriction (Lk)(z) = −1 for z = −a is very weak, because, ifĝ ∈ L 1 (R + , exp a ), then (Lk)(z) is analytic for z > −a and by the Riemann-Lebesgue lemma vanishing as |z| → ∞. Hence we can choose a smaller a such that |a − a | is arbitrarily small and the theorem applies.
The eigenmodes of the Volterra equation can be related to eigenmodes z λ,j of the linear evolution (Equation (5)) on C(N × R). These eigenmodes are vanishing except in the first spatial modes where they satisfy z λ,j (1, ·) ∈ L ∞ (R + , exp a ).
For the center-unstable manifold reduction we therefore look for solutions in sup ξ∈R e aξ |u(l, ξ)|.
The nonlinearity is not controlled in Z a but within By a spectral analysis we can find a continuous projection P cu from Y a to Z a cu := z λi,j : i = 1, . . . , n and j = 0, . . . , p i − 1 with complementary projection P s mapping Z a to Z a s and Y a to Y a s . The image of P s is invariant under the linear evolution and decaying with rate a. In fact the higher modes decay quicker in this norm, so that the solution to the linear evolution with forcing in Y a s is within Z a s . Hence the theory of center-unstable manifold reduction [21,37,38] applies and we can reduce the dynamics for the local behaviour. This is our third main result and here we also consider K as variable in order to discuss the asymptotic result for couplings K close to a fixed coupling K c .
Theorem 7. Given a > 0 and g ∈ M(R) withĝ ∈ L 1 (R + , exp a ). Let K c be a coupling constant such that Theorem 6 applies with center-unstable modes λ 1 , . . . , λ n satisfying cu } is invariant and exponentially attractive under the nonlinear evolution with coupling constant K = K c + in the region {u ∈ Z a : u Z a ≤ δ}. Moreover, the derivatives of ψ can explicitly be computed at 0.
Hence we can reduce the dynamics around the spatially homogeneous distribution to the finite dimensional dynamics of M . We demonstrate the application by recovering Chiba's bifurcation result for the Gaussian distribution [10]. In particular, it shows the nonlinear stability of the partially locked states close to the bifurcation.
The key point for the center-unstable manifold reduction is the use of the regularity in the natural frequency ω. This allows the use of the weighted space Z a on which the linear generator does not have a continuous spectrum along the imaginary axis as seen in earlier linear analysis [35]. Instead the linear evolution decays with rate al in the lth spatial mode except in the first mode l = 1 where a discrete spectrum can arise. Thus the effect of using Z a is very similar to adding white noise [32] of strength D to the system and in this case D > 0 the stability and center manifold reduction have been done [7,15,35].
We remark that the use of Z a is compatible with adding noise and thus for sufficiently regular velocity distribution the center-unstable manifold reduction can be used to study the bifurcation behaviour with the noise strength D as additional parameter in the limit D → 0. The behaviour of the amplitude equations in the center manifold reduction in this limit is reviewed in [34,Section 11] and this explains why taking the limit D → 0 gave the correct behaviour for regular distributions [15]. Note that even in this case the order parameter can decay faster than rate D through the Landau damping mechanism.
In order to understand the norm Z a (which is also the measure for the perturbations in Theorem 4) we relate the norm to the analytic continuation in the strip {z ∈ C : 0 ≤ z ≤ a} of the complex plane. Here we can demonstrate explicitly that for suitable velocity marginals the partially locked states have finite norm in Z a and are small perturbations if their order parameter η is small. In particular, this shows that partially locked states with a small η and coupling K just above K c result in exponential damping when the coupling is lowered below K c by Theorem 4. This is particularly interesting, because these partially locked states are not regular in general, e.g. they do not even have a density.
The exponential stability of the incoherent state and the convergence in Z a to the reduced manifold do not imply convergence in L 2 . Nevertheless, we can relate it to weak convergence for sufficiently nice test functions because the Fourier transform is bounded.
In case of the algebraic stability (Theorem 5), the convergence is very weak but we can conclude ∞ 0 |η(t)|dt < ∞. Going back to the original equation this shows control in the gliding frame and weak convergence of the phase marginal.
The plan of the paper is as follows: In Section 2 we review the notion of the considered solutions and their existence and uniqueness. In Section 3 we show the global energy estimate (Theorem 1). The linear evolution and its stability is discussed in Section 4. In Section 5 we use the linear stability to show nonlinear stability. In preparation for the center-unstable manifold reduction, the linear decomposition is derived in Section 6. With this the nonlinear reduction (Theorem 7) is shown in Section 7. The discussion of the norm and the convergence are in Section 8.
1.3. Acknowledgements. I would like to thank Clément Mouhot for many helpful suggestions and discussions. Furthermore, I would like to thank Bastien Fernandez, David Gérard-Varet and Giambattista Giacomin for inviting me to Paris and discussing a first draft of the work.

Mean-field limit
Early studies of the mean-field limit for the Kuramoto equation are [33,35]. Lancellotti [24] noted that the theory developed for the Vlasov equation [8,16,27] applies because the interaction is Lipschitz. This setup is also reviewed in [9] and was derived using moments in [11].
Since we want to be able to consider partially locked states, we formulate solutions in terms of probability measures. As mentioned in the introduction, we assume the initial data for time t = 0 and solve for t ∈ R + , because we have global existence and uniqueness. Definition 8. Let C M be the solution space C w (R + , M(Γ)), which consists of the families of weakly continuous probability measures on Γ, i.e. ρ · ∈ C M is a family The trajectories T t,s [ρ · ] : Γ → Γ of the system are defined as T t,s [ρ · ](Q) = P (t) where P : R + → Γ is the solution to the initial value problem d dt P (t) = (V [ρ · ](t, P (t)), 0), P (s) = Q.
Since |η(t)| ≤ 1 for all t ∈ R + , the velocity field V [ρ · ] is globally Lipschitz and the trajectories are well-defined. Moreover, all derivatives of V [ρ · ](t, P ) with respect to P are bounded, so that T t,0 [ρ](P ) also has uniformly bounded derivatives of all orders for finite ranges of t. Finally, T t,s is invertible with inverse T s,t and thus a diffeomorphism.
By standard arguments for scalar transport equation, the following lemma holds. The well-posedness implies that this is the mean-field limit as N → ∞. For this consider the empirical measure ρ N,t of a system of N oscillators and suppose ρ N,0 ρ in . As the empirical measure is a solution, the well-posedness shows that at later times t also ρ N,t ρ t . Lemma 9 shows that the regularity is preserved along the flow as we can see from the following lemma.
Lemma 11. Let ρ in ∈ M(Γ) and ρ · ∈ C M be the solution to the Kuramoto equation with initial data ρ in . If ρ in has a density f in , then for every t ∈ R + the measure ρ t has a density f (t, ·) given by the flow Lemma 9 also shows that the velocity marginal is constant in time.
Lemma 12. Let ρ in ∈ M(Γ) and g ∈ M(R) be its velocity marginal. If ρ · ∈ C M is the solution to the Kuramoto equation, then g is for all t ∈ R + the velocity marginal of ρ t .
For a solution ρ · ∈ C M and velocity distribution g ∈ M(R) define the Fourier transform u respectivelyĝ by Equations (1) and (2). The evolution of the difference to the uniform state is captured by the restriction to l ≥ 1. Its evolution is described in Fourier variables as follows: Theorem 13. Suppose ρ in ∈ M(Γ) with velocity marginal g ∈ M(R). Let ρ · ∈ C M be the solution to the Kuramoto equation. Then the Fourier transform u is in Proof. For every t ∈ R + , the distribution ρ t is a probability measure, so that u(t, l, ξ) is bounded by 1 and continuous with respect to ξ. Since ρ · is weakly continuous, u is also continuous with respect to time.
By Lemma 12 the velocity marginal is always g so that Plancherel's formula shows that Equation (11) The free transport is a well-posed linear problem, so that Duhamel's principle on the free transport holds.
Proof. In the weak formulation Equation (11) consider all terms except the free transport as forcing. Then the problem has a unique solution given by v.
For the stability and bifurcation result we construct solutions u to Equation (11). By the following theorem these correspond to the Fourier transform of the meanfield solution. Then by Lemma 14 for t ∈ [0, T ], l ∈ N and ξ ∈ R we have e(s).
Therefore, we must have e(t) = 0 for t ≤ t * , i.e. the solutions agree. Since we can repeat the argument in steps of t * up to time T , this shows uniqueness up to time T .
By Theorem 13 the Fourier transform u of a solution to the Kuramoto equation satisfies the assumptions of the previous uniqueness theorem for every α and β.
Theorem 16. The restriction of the Fourier transform u to C(R + × N × R + ), i.e. ξ ≥ 0, satisfies an appropriate weak PDE (Theorem 13) to which the Duhamel formula (Lemma 14) applies and the uniqueness holds for the restriction ξ ≥ 0 in Theorem 15.
Proof. Characteristics of the full equation are never entering the restricted region and the nonlinear interaction in the restricted region is determined by the region itself. Hence the restriction satisfies the appropriate restriction of the weak PDE Theorem 13 whose solution is given by the Duhamel formula Lemma 14. The proof of the uniqueness theorem holds as before when restricting ξ to R + everywhere.

Global stability by energy method
The basic idea is to differentiate the energy functional in time and to note that the nonlinear terms cancel, because they are skew-Hermitian. However, the Fourier transform u does not need to be differentiable. Therefore, we use an approximation procedure.
For the Fourier transform u ∈ C(R + × N × R) of a solution to the Kuramoto equation define the functional I(t) by Equation (4). Then for all t ∈ R + Proof. In order to overcome the non-differentiability, let χ ∈ C ∞ c (R) be a nonnegative function with x∈R χ(x)dx = 1 and let χ δ (x) = δ −1 χ(x/δ) for δ > 0. Define the mollified u δ by Since u is bounded by 1, also u δ is bounded by 1 and has a bounded derivative with respect to ξ. Moreover, u δ (t, l, ·) 2 ≤ u(t, l, ·) 2 .
By Lemma 14 whereĝ δ and u in,δ are the mollifications ofĝ and u in , respectively. Hence u δ is also continuously differentiable with respect to t ∈ R + and satisfies classically and for l ≥ 2 For , ζ > 0 define Since u δ and ∂ ξ u δ are bounded, I δ, ,ζ (t) is differentiable and we can differentiate under the integral sign to find As |η(t)| ≤ 1, this shows by Gronwall's inequality for a fixed time t Taking δ → 0, we find by dominated convergence The integral can be controlled as Hence, the previous bound shows The global stability result, Theorem 1, follows from this lemma by finding a suitable weight φ such that α < 1. For this we want to minimise the integral defining α. Formally, the Euler-Lagrange equation is With φ(0) = 1 the solution formally is ∞ 0 |ĝ(ξ)|dξ, this indeed defines a bounded increasing weight and we find Thus if K < K ec , for a suitable A we have the required control.
Thus φ is a valid weight for the theorem and by Lemma 17 a solution satisfies the required bound For the initial bound, use Plancherel to note

Stability of the linearised evolution
As in the mean-field limit, the linear evolution given by Equation (5) is understood as PDE for continuous solutions in the same weak sense, i.e. tested against C 1 0 functions. Proof of Lemma 2. The free transport has a unique weak solution, so that, as in Lemma 14, Duhamel's formula implies and for l ≥ 2 (e tL u in )(l, ξ) = u in (l, ξ + lt). Taking ξ = 0 in Equation (12) shows that ν(t) = (e tL u in )(1, 0) must satisfy Equation (6). By the general theory of Volterra equation [19, Chapter 2] Equation (6) has a unique solution and the solution is continuous. Hence the evolution operator e tL is well-defined and given by Equations (6), (8) and (9).
The solution of the Volterra equation can be characterised by the following lemma.
Lemma 18. For a convolution kernel k ∈ L 1 loc (R + ) there exists a unique resolvent kernel r ∈ L 1 loc (R + ) satisfying r + r * k = k. For f ∈ L 1 loc (R + ) the unique solution x ∈ L 1 loc (R + ) to the Volterra equation x + k * x = f is given by x = f − r * f . Proof. This is Theorem 3.1 and 3.5 of Chapter 2 of [19].
We can control the convolution by using the (sub)multiplicative structure of the weights.
and for A ≥ 1 and b ≥ 0 showing the first bound.
For A ≥ 1 and b ≥ 0 and t ∈ R + Taking the supremum shows the second claim.
Proof. By the previous lemma and the explicit formula for ν, Equation (8) shows for ξ ∈ R |(e tL u in )(1, ξ)|e a(t+ξ) ≤ u in L ∞ (R + ,exp a ) which proves the claim.
The stability condition on the kernel comes from the literature.
Proof of Theorem 3. The exponential condition follows from Lemma 3.4 and Theorem 4.1 of Section 2 of [19] and is called the Paley-Wiener theorem. The algebraic condition comes from Gel'fand's theorem as Theorem 4.3 in Section 4 in [19], because the weight p b is submultiplicative.
The stability criterion 1 + Lk(z) = 0 for z ≥ 0 can be understood as imposing that Lĝ({z : z ≥ 0}) does not contain 2/K. Since Lĝ is analytic, the argument principle shows that this is equivalent to imposing that Lĝ(iR) does not encircle 2/K.
For z > 0 we find by Fubini By continuity of Lĝ we thus find where the last equality is Plemelj formula for continuous g and PV denotes the principal value integral. The curve Lĝ(iR) starts and ends at 0, is bounded, and goes through the right half plane where it crosses the real axis by continuity of the principal value integral. If g > 0, this shows that it crosses the real axis at a positive value. Thus there exists a critical value K c such that for 0 ≤ K < K c the system in linearly stable and for K c + δ > K > K c with some δ > 0 the system is unstable.
For the Gaussian distribution, we can plot the boundary explicitly (cf. Figure 1) and note that for K < 2/(πg(0)) the solution is stable while for K > 2/(πg(0)) there is exactly one root of (1 + Lk)(z) for z ≥ 0.  Already for two Gaussians with variance 1 centred at ω = ±1.5, we see a more interesting behaviour (cf. Figure 2). Here we see that at the critical coupling two roots of (1 + Lk)(z) appear and that for sufficiently large K this reduces to one root.
A sufficient condition for stability is In this case (Lĝ(ix)) < 2/K so that Lĝ(iR) cannot encircle 2/K and the solution is stable. In the case of a symmetric distribution with a maximum at 0, this condition is sharp for the first instability because then (Lĝ)(0) = πg(0). Comparing the linear stability condition to possible stationary solutions, we can understand the stability condition on the real part πg(x) as critical mass density to form an instability while the imaginary part ensures that it is not drifted away.

Nonlinear stability
In Fourier variables the nonlinear interaction is given by R : C(N×R) → C(N×R) defined as Here we see that the nonlinearity is not bounded as it increases with the spatial mode l. In the exponential convergence we can match this with the faster decay of the linear operator in higher modes. In the algebraic case we compensate for this in the norm by the factor (1 + t) −α(l−1) .
We characterise the solution by Duhamel's principle on the linear evolution.
Proof. For u consider the nonlinear terms as forcing in Equation (11). The linear problem has a unique solution given by T u. Hence u is a solution if and only if T u = u.
For the exponential and algebraic stability we want to propagate the norms given in the introduction using Lemma 21. However, we do not know a priori if these norms stay finite. Therefore, we construct a solution using Banach fixed point theorem and by the uniqueness (Theorem 15) this must be the mean-field solution.
By the assumption Lemma 20 applies with a finite constant c 1 . This shows for u ∈ C(R + × N × R) and t ∈ R + , ξ ∈ R If u,ũ ∈ C(R + × N × R) satisfy u e , ũ e ≤ (1 + c 1 )δ, then as before we find Hence we can choose δ small enough such that T is a contraction on The structure for the proof of algebraic stability is very similar to the previous proof of exponential estimate.
Proof of Theorem 5. Define on C(R + × N × R) the norm · a by for u ∈ C(R + × N × R). By the assumption Lemma 20 applies with a finite constant γ 1 . This shows for For modes l ≥ 2 the explicit formula Lemma 2 shows For l ∈ N with αl > b − 1 we can bound by a finite constant independent of l.
For l ∈ N with l ≥ 2 and αl − b ≤ 1 the integral t 0 (1 + s) αl−b ds grows at most logarithmically, so that it can be bounded by (1 + t) α(l−1) . Hence there exists a finite γ 2 satisfying For two functions u,ũ ∈ C(R + × N × R) with u a , ũ a ≤ (1 + γ 1 )δ, we find for t ∈ R + and ξ ∈ R + as before and for l ≥ 2 The decomposition now follows from Theorem 2.1 of Chapter 7 of [19].
For the characterisation of the stable subspace we use the functional α.
Lemma 22. For a ∈ R and λ ∈ C and j ∈ N with λ > −a define the continuous functional α λ,j on Y a by Proof. Since |u(1, ξ)| ≤ e −aξ u Y a , the functional is well-defined, continuous and for j ≥ 1 the derivative can be computed by differentiating under the integral sign.
Conversely, assume ν ∈ L ∞ (R + , exp a ). Then at z ∈ C with z > −a the Laplace transform of ν, k and u(1, ·) are analytic and by the Volterra equation (6)  Hence the LHS has roots λ 1 , . . . , λ n of multiplicities p 1 , . . . , p n . Thus the RHS must have the same roots, i.e. α λi,j (u) = 0 for i = 1, . . . , n and j = 0, . . . , p i − 1, which shows u ∈ Y a s . The growth condition ν ∈ L ∞ (R + , exp a ) under the linear evolution is invariant, so that Y a s is invariant under the linear evolution. We can use this to control the decay in Z a s = Y a s ∩ Z a under the linear evolution and the forcing term from a Duhamel equation.
Lemma 24. Assume the hypothesis of Lemma 23 with a ≥ 0. Then for u ∈ Z a s we Remark 25. The forcing F is measured in the weaker norm Y a and we are able to recover the control in the stronger norm Z a .
Proof. Since u is in Z a s , we find by Lemma 23 for t ∈ R + and ξ ∈ R e aξ |(e tL u)(1, ξ)| ≤ c 1 e −at u Z a and by (9) for l ≥ 2 e aξ |(e tL u)(l, ξ)| ≤ e −alt u Z a which shows the first claim.
These bounds are uniformly integrable as Hence the integral is well-defined and defines a continuous function v ∈ C(N × R). Moreover, it shows the claimed control of v Z a .
In order to find the eigenmodes in Z a and the projection to them, we consider the linear generator L on Z a . Lemma 26. Assume the hypothesis of Lemma 23 with a > 0. Define the closed operator L : and for l ≥ 2 by (Lu)(1, ξ) = l∂ ξ u(l, ξ).
Since λ ∈ σ, this determines u(1, 0) uniquely and the solution is given by Equation (14). Moreover, the found solution satisfies u Z a < ∞. Hence such a λ is not in the spectrum of L and the resolvent map takes the given form.
If λ ∈ σ, the above shows that L is not injective, so that σ is in the spectrum.
Since we assumeĝ ∈ L 1 (R, exp a ), we always have z λ,j ∈ D(L) ⊂ Z a .
Using holomorphic functional calculus, we can construct the required projection.
Lemma 28. Assume the hypothesis of Lemma 23 with a > 0. There exists a continuous projection P cu : Y a → Z a cu such that the complementary projection P s = I − P cu maps to Y a s . Proof. Take γ a contour in {z ∈ C : z > −a} encircling all roots λ 1 , . . . , λ n once and define Then P cu is a continuous projection. Since its value is given by the residues, Equation (14) shows that P cu maps into Z a cu . On the image of the complementary projection P s the resolvent map (L − λ) −1 is continuous for λ > −a. By Equation (14) this can only be true if for all roots λ 1 , . . . , λ n with multiplicities p 1 , . . . , p n also Lu(1, ·) has a root of matching multiplicity, i.e. α λi,j must be vanishing for i = 1, . . . , n and j = 0, . . . , p i − 1. Hence the image is within Y a s .

Center-unstable manifold reduction
Since the nonlinearity R is smooth from Z a to Y a , the center-unstable manifold reduction follows from the decay on the stable part (Lemma 24) and the continuous projections (Lemma 28). The center-manifold reduction in Banach spaces is discussed in [38,Chapter 2], which can be adapted to our case with a weak PDE as Duhamel's principle holds again. Moreover, instead of a center manifold reduction we do a center-unstable manifold reduction [37, Section 1.5] as this implies convergence to the reduced manifold. The statements are also discussed in [21].
In order to understand the bifurcation behaviour, the center-unstable manifold is constructed with = K − K c as additional parameter. Furthermore, we need to localise the non-linearity around the incoherent state. This localisation might affect the resulting reduced manifold, which is not unique. Moreover, we choose the localisation to preserve the rotation symmetry of the problem so that the resulting manifold has the same symmetry.
Hence we study the solution inZ a = Z a × R andỸ a = Y a × R with norms (u, ) Za = max( u Z a , | |) and (u, ) Ỹa = max( u Y a , | |). In these define the linear evolution operator e tL as e tLũ = e tL (u, ) = (e tL u, ).
The center-unstable spaceZ a cu is now spanned by (0, 1) and (z λi,j , 0) for i = 1, . . . , n and j = 0, . . . , p i − 1. The projectionP cu is defined asP cu (u, ) = (P cu u, ) with P cu from Lemma 28. The complementary projectionP s maps into the stable part Y a s = Y a s × {0}. For the center-unstable manifold reduction the used properties of the linear evolution are collected in the following lemma. the linear evolution acts as finite-dimensional matrix exponential whose generator has spectrum λ 1 , . . . , λ n . OnZ a s =P sZ a there exists a constant c 1 such that for t ∈ R + and u ∈Z a s e tL u Za ≤ c 1 e −at u Za and there exists a continuous function c : [0, a) → R + such that for µ ∈ [0, a) and Proof. The linear evolution on is constant so that the properties follow from the study of e tL , i.e. Lemmas 24, 27 and 28.
We now study the evolution PDE forũ = (u, ) ∈ C(∆ × N × R × R) understood weakly against test functions in C 1 0 and where ∆ ⊂ R is a time interval. (16) satisfies Duhamel's formula, i.e. for t 0 , t 1 ∈ ∆ with t 0 ≤ t 1 the following holds

Lemma 30. A solutionũ of Equation
Proof. ConsiderR δ (ũ) as forcing. Then the linear problem has a unique solution given by the formula.
Lemma 31. Letũ be a solution of Equation (16). Givenũ(t 0 ) ∈Z a the solution is uniquely determined at later times t ≥ t 0 and ũ(t) Za is uniformly bounded for compact intervals of later times. If δ is small enough, there exists, givenũ(t 0 ) ∈Z a , a global solution for t ≥ t 0 .
Proof. The projections commute with the linear evolution in the Duhamel formula and are bounded. SinceR δ is bounded, this shows with the control of the linear evolution (Lemma 29) that ũ(t) Za is uniformly bounded for compact time-intervals of later times. Furthermore,R δ is Lipschitz so that the Duhamel formulation gives a Gronwall control showing uniqueness.
If δ is small enough, then the Picard iteratioñ is a contraction, so that a solution exists.
With this we can construct the reduced manifold. We want to capture the evolution in the reduced manifold where the stable part has decayed, which we can characterise by the growth as t → −∞.
The growth condition already suggest that this captures the dominant behaviour and M is attractive.
Proof. Since we have a Duhamel formula, the proofs of Theorem 1 and 2 of [38] work. We further adapt it, as in Section 1.5 of [37], to the center-unstable case. Parts of the results are also stated as Theorem 2.9 and 3.22 of Section 2 of [21].
The key point is to note that by the Duhamel formulation a solution with the growth bound is a fixed point of the map T given by for t ∈ R − on the space C(R − ×N×R×R) with norm ũ = sup t∈R − e µt ũ(t, ·) Za . Here note that onZ a cu the linear evolution e tL restricts to a finite dimensional matrix exponential, which has a well-defined inverse e −tL .
For small enough δ the map T is a contraction with a unique fixed point depending onP cuũ0 so that there existsψ :Z a cu →Z a s such that the unique solution isP cuũ0 +ψ(P cuũ0 ) [38,Theorem 1]. Hence the set M is a reduced manifold with the given form. Moreover, the growth condition is invariant under the evolution.
By a Fibre contraction argument ( [38, Theorem 2] and [37, Section 1.3]), for small enough δ, the mapψ has the claimed regularity and the derivatives can be computed explicitly at 0 with the given form.
As discussed in [37, Section 1.5] and [21, Theorem 3.22 of Chapter 2] a similar fixed point argument shows the exponential attractiveness.
Hence, with a sufficient localisation, we can prove Theorem 7.
Proof of Theorem 7. Apply Theorem 32 with a small enough δ. Then for | | ≤ δ and u ≤ δ the nonlinearity ofR δ agrees with the nonlinearity of the Kuramoto equation, i.e. a solution to the localised evolution is also a solution to the original problem in the restricted region.
The mapψ maps intoZ a s and thus has the formψ(u, ) = (ψ(u, ), 0) for ψ mapping into Z a s . As is constant, the reduced manifold must have the claimed form and the result follows from the the uniqueness (Theorem 15).
Remark 33. The Kuramoto equation has the rotation symmetry, i.e. if u is a solution to the Kuramoto equation, thenū(t, l, ξ) = e ilα u(t, l, ξ) is again a solution. This symmetry is also satisfied by the localised nonlinearity and thus by the reduced manifold.
7.1. Bifurcation analysis. The center-manifold reduction can be used to determine the bifurcation behaviour by studying the evolution on the reduced manifold. The resulting behaviour crucially depends on the velocity distribution and [15] contains several examples based on the center-manifold reduction with noise which at this point is very similar. The Penrose diagrams (cf. Figures 1 and 2) already show the dimension of the reduced manifold as the covering number, e.g. this shows for the example of Figure 2 that at the critical coupling two eigenmodes appear.
As an example we repeat Chiba's analysis [10] for the bifurcation of the Gaussian distribution. In this example the density is given by and the Fourier transform isĝ From the discussion in Section 4, we see that the critical coupling is and that at the critical coupling K c there exists a single eigenvalue at λ = 0 for small enough a.
Recall, that the eigenvector z 0,0 is given by Hence z 0,0 (1, 0) = 1 and Now apply Theorem 7 with sufficiently small a > 0 such that λ 1 = 0 is the only eigenmode and with k ≥ 3. Then forũ = (u, ) ∈Z a cu ψ(ũ) For this we find that only the first two spatial modes are non-vanishing with Hence we can compute Hence we find On the reduced manifold the solution is β(t)z 0,0 + ψ(β(t)z 0,0 , ), which evolves by the previous computation as For small enough ≤ 0, thus β = 0 is a stable fixed point, i.e. the incoherent state is stable. For small enough > 0, the zero solution is unstable, but the set |β| = β c are stable attractors where Hence we have proved the claimed bifurcation from the incoherent state and that for small enough the appearing stationary states are nonlinearly stable. Since z 0,0 (1, 0) = 1, the value of β agrees with the order parameter so that we have found the known result.

Boundedness and convergence in the exponential norms
The exponential norm L ∞ (R, exp a ) can be related to an analytic extension in a strip {z : 0 ≤ z < a}. exp a ). Then f has an analytic continuation to {z : 0 ≤ z < a}.
Proof. For z ∈ R the Fourier inversion formula shows By the assumed supremum bound onf for ξ ≥ 0 and the bound |f (ξ)| ≤ f 1 for ξ ≤ 0, this definition extends to the given range and has a complex derivative given by differentiating under the integral sign.
Lemma 35. Let f ∈ L 1 (R) and suppose f has an analytic continuation to {z : 0 ≤ z ≤ a} such that |f (z)| → 0 uniformly as | z| → ∞ and f (· + ia) ∈ L 1 (R). Then the Fourier transformf is bounded as Proof. By the assumed decay we can deform the contour of integration aŝ Bounding the integral with the L 1 norm of f shows the claimed result.
We can use this to show the finite Z a norm of partially locked states (see [34,Chapter 4] for a brief review of these states).
Theorem 36. Suppose a velocity marginal with density g ∈ L 1 (R) and that g has an analytic continuation to {z : 0 ≤ z ≤ a} such that |g(z)| → 0 uniformly as | z| → ∞ and g(· + ia) ∈ L 1 (R). Then a partially locked state ρ locked has Fourier transform u ∈ Z a . If furthermore |η| ≤ a/( √ 2K), then The assumption on the velocity distribution is for example satisfied for a Gaussian distribution. This estimate shows that partially locked states with small order parameter η are small perturbations in Z a .
Proof. Suppose a partially locked state with order parameter η. Then for |ω| < K|η| all oscillators are trapped at θ + satisfying For |ω| > K|η| the oscillators are distributed proportionally to the density .
By Cauchy's residue theorem f (l, ω) = g(ω) η |η| β ω K|η| l where β is given by Now β has an analytic continuation to the complex upper half plane and is bounded there by 1. Hence Lemma 35 shows for all l ∈ N u(l, ·) L ∞ (R,exp a ) ≤ g(· + ia) L 1 (R) , which implies that u ∈ Z a . For |z| ≥ √ 2 we can estimate by the series expansion Hence we can estimate for |η| ≤ a/( √ 2K) that which is the claimed bound.
We can relate convergence in Z a as found in the center-unstable manifold reduction and the exponential stability with weak convergence.
Theorem 37. Let ρ · ∈ C M be a solution to the Kuramoto equation with initial data ρ in ∈ M(Γ) and Fourier transform u. If u(t, ·, ·) → v in Z a as t → ∞ where v is the Fourier transform of ρ ∞ ∈ M(Γ), then for every φ ∈ H 4 (Γ) the integral ρ t (φ) converges to ρ ∞ (φ) as t → ∞.
In the stability theorem 5 we only control a norm in half of the Fourier coefficients which also looses control for higher spatial modes. However, this controls the order parameter η and in particular it shows ∞ 0 |η(t)|dt < ∞. With this knowledge we can go back to the original equation to deduce convergence properties. As an example we show a convergence in the gliding frame.
In the gliding frame the position of each oscillator is corrected by the effect of the free transport, i.e. if ρ · is the solution of the Kuramoto equation which has a density f , then the density h in the gliding frame is given by h(t, θ, ω) = f (t, θ + tω, ω).
If f in ∈ L 2 (Γ) and ∂ θ f in ∈ L 2 (Γ), we can integrate the previous inequality over ω to show the convergence of h(t, ·, ·) in L 2 (Γ) as t → ∞.
The convergence in the gliding frame implies for example weak convergence in the normal setup.
Given > 0, we then have for large enough t by the convergence in L 2