Uniform stability and mean-field limit for the augmented Kuramoto model

We present two uniform estimates on stability and mean-field limit for the "augmented Kuramoto model (AKM)" arising from the second-order lifting of the first-order Kuramoto model (KM) for synchronization. In particular, we address three issues such as synchronization estimate, uniform stability and mean-field limit which are valid uniformly in time for the AKM. The derived mean-field equation for the AKM corresponds to the dissipative Vlasov-McKean type equation. The kinetic Kuramoto equation for distributed natural frequencies is not compatible with the frequency variance functional approach for the complete synchronization. In contrast, the kinetic equation for the AKM has a similar structural similarity with the kinetic Cucker-Smale equation which admits the Lyapunov functional approach for the variance. We present sufficient frameworks leading to the uniform stability and mean-field limit for the AKM.


1.
Introduction. Synchronization of weakly coupled oscillators is ubiquitous in our nature, e.g., rhythmic heart beatings of pacemaker cells, synchronous flashing of fireflies and collective hand clapping in a concert hall, etc [1,6,34,35]. After Huygen's observation on two pendulum clocks hanging on the same bar, collective behaviors of weakly coupled oscillators have been reported from time to time in scientific literature (see [34]). However, major scientific progress on the collective dynamics of complex systems was initiated by Winfree and Kuramoto about a half century ago in [26,27,41]. Recently, research on the collective dynamics of complex systems has received lots of attention due to engineering applications in sensor network, mobile network and control of unmanned aerial vehicles (UAV) etc. In [22], the authors observed a formal analogy between the Cucker-Smale flocking model and the Kuramoto model for synchronization, and provide a quantitive estimate for the synchronization based on Lyapunov functional approach. In the paper, we further investigate this formal analogy and study dynamic asymptotic properties of the AKM.
Consider an ensemble of Kuramoto oscillators lying on the nodes of the complete graph with N -nodes, and assume that the state of an oscillator is described by a real-valued function "phase". Let θ i be the phase of the i-th oscillator whose dynamics is given by the Kuramoto model: where κ and ν i denote the coupling strength, and the natural frequency of the i-th Kuramoto oscillator, respectively. On the other hand, it is well known in [28] that the dynamics of system (1) with N 1 can be described by the corresponding mean-field equation, namely kinetic Kuramoto equation. To be specific, let f = f (θ, ν, t) be the one-particle distribution function at phase θ, natural frequency ν at time t. Then, the kinetic Kuramoto equation reads as follows.
Next, we briefly explain how to bypass the aforementioned difficulty for the complete synchronization of the kinetic equation (2) for distributed natural frequencies.
Our idea is to lift the first-order model (1) into a second-order model by introducing an auxiliary frequency variable ω i =θ i , i.e., we differentiate the equation (1) with respect to t and obtain an augmented Kuramoto model (AKM): On the other hand, we consider a mean-field limit N → ∞. In this case, we set f = f (θ, ω, t) to be a probability density function corresponding to (4). Then, a formal BBGKY Hierarchy argument yields the kinetic equation: As an analogy with the kinetic Cucker-Smale model in [16], we can use the frequency variance functional to measure the emergence of complete synchronization for (5).
In this paper, we are interested in the following questions for (4) and (5): Under what conditions, can system (4) exhibit the complete synchronization? • (Q2): Is the system (4) uniformly p -stable with respect to initial data? • (Q3): Can we derive the mean-field kinetic equation (5) from the particle model (4) as N → ∞ uniformly in time?
The first two questions might be generalized to the locally coupled Kuramoto model on a general symmetric and connected networks. However, the last question, i.e., uniform mean-field limit can be treated only for mean-field couplings (e.g., BBGKY hierarchy arguments break down for the locally coupled case). As aforementioned, since our main motivation is to study the complete synchronization of the kinetic level in a direct manner, we consider only the complete network case.
The main results of this paper are four-fold: First, we provide a sufficient framework for the complete synchronization estimate. Our sufficient conditions are expressed in terms of the coupling strength κ, the diameter of the set of natural frequencies and initial data, and they are free of the number of oscillators (see Theorem 3.2). Second, we provide the uniform stability estimate of (4) with respect to initial data in a metric equivalent to the p -distance in phase space. Our uniform stability roughly says that the p -distance between two configurations at time t is uniformly bounded by the constant multiple of initial p -distance between two initial data (see Theorem 4.3). Third, we present a uniform-in-time mean-field limit for (4) as a direct application of the exponential flocking estimate in Theorem 3.2 and uniform-in-time stability estimates in Theorem 4.3. Our last result is to derive the complete synchronization estimate for the mean-field kinetic equation (5) using a robust Lyapunov functional approach.
The rest of this paper is organized as follows. In Section 2, we briefly discuss theoretical minimum for the Kuramoto model and its augmented model. In Section 3, we present a synchronization estimate for the AKM (4). In Section 4, we present a uniform p -stability estimate for the augmented model. In Section 5, we study the uniform mean-field limit from the particle model (4) to the corresponding kinetic equation uniformly in time, and we also study the complete frequency synchronization estimate for the kinetic equation. In Section 6, we provide a synchronization estimate for the kinetic Kuramoto equation (5) with distributed natural frequencies.
Finally, Section 7 is devoted to a brief summary of our main results and discussion for future works.
Before we proceed to the next section, we introduce the notation which will be used in the rest of the paper.
In the following discussion, we only consider the case in which the oscillators are confined in a half circle. It is obvious that |θ| o = |θ| for θ ∈ (−π, π).

Preliminaries.
In this section, we briefly review a theoretical minimum for the Kuramoto model and the augmented Kuramoto model.
2.1. The Kuramoto model. In this subsection, we briefly discuss an associated conservation law, and review the state-of-the-art results on the complete synchronization for the Kuramoto model. First, we introduce a time-dependent quantity C(Θ, V, t):

301
Lemma 2.1. Let Θ = Θ(t) be a phase vector whose dynamics is governed by (1). Then, the quantity C(Θ, V, t) is a constant of motion.
Proof. Let Θ = Θ(t) be a Kuramoto flow. Then, we have This yields the desired estimate.
Next, we discuss the equilibrium for the Kuramoto model (1). Note that the equilibrium solution Θ = (θ 1 , · · · , θ N ) is a solution to the following equilibrium system: We sum up the equation (6) with respect to i to obtain Thus, if the system (6) does have a solution, then the total sum of natural frequencies is zero. Hence if N i=1 ν i = 0, then system (6) does not have a solution. This leads to the need of the relaxed equilibria in the following definition. We first recall several definitions for a phase-locked state and asymptotic phase-locking.
1. Θ is a phase-locked state if all relative phase differences are constants: 2. Θ exhibits asymptotic phase-locking (complete synchronization) if the relative frequencies tend to zero asymptotically: Remark 2. Note that solutions to the following equilibrium system: correspond to phase-locked states of (1).
We next briefly review the state-of-the-art result for the Kuramoto model (1). It is well known in [11,24,36] that system (1) can be lifted as a dynamical system on R N and can also be written as a gradient flow with an analytical potential in R N : For a gradient flow system with analytical potential, uniform boundedness is equivalent to the convergence of solution toward the phase-locked state. As a dynamical system in R N , the uniform boundedness of (1) in R N in a rotating frame moving with the speed of average of natural frequencies is not obvious since nonidentical oscillators can cross each other. In the following theorem, we summarize the state-of-the-art results for the emergence of phase-locked states from generic initial data in a large coupling strength regime with κ 1.
1. Suppose that the initial phase configuration Θ 0 is confined in a half circle and the coupling strength κ is positive such that Then, for any solution Θ = Θ(t) to (1), we have lim t→∞ D(Θ(t)) = 0.
2. Suppose that the initial phase configuration Θ 0 is confined in a half circle and the coupling strength κ is sufficient large such that Then, for any solution Θ = Θ(t) to (1), there exists a positive constant λ such that 3. Suppose that natural frequencies are distributed and initial configurations satisfy Then there exists a large coupling strength κ ∞ > 0 such that if κ ≥ κ ∞ there exists a phase-locked state Θ ∞ such that the solution with initial data Θ 0 satisfies Remark 3. 1. In the course of the proof for the first statement, we can show that there exists a finite time t 0 and D ∞ ∈ (0, π 2 ) such that 2. The result of [11] does not yield detailed asymptotic dynamics of identical Kuramoto oscillators. However, when the diameter of the emergent phase-locked state is less than π and the coupling strength is sufficiently large, then we can show that convergence speed is at least exponential. See [9,12,13] for detailed discussion.
2.2. The augmented Kuramoto model. In this subsection, we discuss basic structural properties of the AKM and relationship between the Kuramoto model and AKM. We set averaged quantities and fluctuations of phase and frequency around them: Then, it is easy to see that the averaged quantities and fluctuations satisfẏ be a solution to (4). Then, the averaged quantities (θ c , ω c ) satisfy the following relations: Proof. We sum (4) with respect to i and use the skew symmetry of cos( The second relation follows from (7) 1 .
Next, we discuss the relation between the first-order model (1) and the secondorder model (4) which is stated in the following theorem. 1 Proof. (i) Let Θ = Θ(t) be a solution to (1) with initial data Θ 0 . Then, it satisfieṡ We set ω i =θ i (9) and differentiate the above equation to obtaiṅ We use (8) to find Letting t → 0+ in (11), we obtain Finally, we combine (9), (10) and (12) to see that Then, we use the relations: to integrate (10) to obtaiṅ Then, we set Finally, we combine (13) and (14) to recover the Kuramoto model.
Before we close this section, we quote the following lemma to be crucially used in later sections.
Lemma 2.6. [21] Suppose that two nonnegative Lipschitz functions X and V satisfy the system of differential inequalities: where α and γ are positive constants. Then, X and V satisfy the uniform bound and decay estimates: where M is a positive constant defined by, 3. Emergence of the complete synchronization. In this section, we present complete synchronization estimates for the AKM (4) in ∞ and p , p ∈ [1, ∞) frameworks by deriving a system of Grönwall's inequalities.
3.1. ∞ -framework. In this subsection, we present a synchronization estimate in ∞ -framework. In each case, the synchronization estimates are obtained in the following three steps: • Step A (Existence of positively invariant set). We identify a positively invariant set which is translation invariant in phase space. • Step B (Derivation of Grönwall's inequality). We introduce a Lyapunov type functional and derive a Grönwall type differential inequality. • Step C (Complete synchronization estimate). Once we derive a Grönwall type inequality for a suitable Lyapunov functional, suitable Grönwall's lemma and continuity arguments yield the desired synchronization estimate.
As candidates for Lyapunov functionals and an invariant set, we introduce phase and frequency diameters: Lemma 3.1. Suppose that initial data (Θ 0 , Ω 0 ) and the coupling strength κ satisfy Then, the set S is positively invariant under the flow (4), i.e., for any solution Θ = Θ(t) with initial data Θ 0 ∈ S, we have Proof. Let Θ 0 be a initial data with D(Θ 0 ) < D ∞ . Suppose that S is not positively invariant under flow. Then, there exists a finite time t * such that By the continuity, we have On the other hand, for any indices i and j, we integrate (4) 2 from 0 to t for t < t * to obtain  Then, for t ∈ [0, t * ], we get For 0 ≤ s ≤ t * , we have Therefore, it follows from (15) and (16) that we have This yields We set Then, it is easy to see thaṫ u(t) = D(Θ(t)), u(0) = 0,u(0) = D(Θ(0)).
Then, the relation (18) is equivalent tȯ Then, (19) and (20) yield On the other hand, since D(Θ(t * )) = D ∞ , there exist indices i and j such that Then, it follows from (4) 1 that we have where we used the hypothesis on κ and (21), which yields contradiction.

UNIFORM STABILITY AND MEAN-FIELD LIMIT OF THE KURAMOTO MODEL 307
Theorem 3.2. Suppose that initial data and coupling strength satisfy Then, we have an exponential synchronization: Proof. Due to the conservation law in Lemma 2.4, we have We set extremal indices M and m such that Then, it follows from (4) 2 that we havė Similarly, we haveω Then, it follows from (22) and (23) This yields the desired exponential decay estimate.
3.2. p -framework with p ∈ [1, ∞). In this subsection, we present p -estimate for (4) for later use. For phase and frequency vectors Θ = (θ 1 , · · · , θ N ) and Ω = (ω 1 , · · · , ω N ), we set Θ p and Ω p : Proposition 1. Suppose that initial data and coupling strength satisfy Then for any solution Proof. (i) Note that We multiply by p|θ i | p−1 to the above relation, take a sum the resulting relation, and use Hölder's inequality to get the following estimate: This yields the desired first differential inequality.
(ii) It follows from (4) 2 that we have We use (25) to obtain We use the monotonicity of f (x) = |x| p−2 x to see Then, we use (26), (27), i=1 ω i = 0 and a priori condition: This yields the desired second differential inequality.
Finally, we combine Proposition 1 and Lemma 3.1 to derive the exponential synchronization. (4) with initial data and coupling strength: Then, there exists a positive constant θ ∞ p such that for p ∈ [1, ∞), The exponential decay of Ω follows from the second equation of (24). On the other hand, it follows from (24) 1 that we have Thus, we have Thanks to Theorem 3.3, we can conclude that there exists a unique phase-locked state Θ ∞ . Moreover, Θ(t) will tend to Θ ∞ exponentially.
Then, for any solution {(θ i , ω i )}, there exists a unique phase lock state Proof. Let Θ = Θ(t) be a solution to system (4). Then, since κ is sufficiently large, we have sup Then, we use Theorem 3.2 to obtain Then for any ε > 0, we can find a positive number T such that ift ≥ T and t ≥ T , then This immediately implies that there exists a unique asymptotic limit θ ∞ i . Moreover, we combine (28) to show that 4. Uniform p -stability estimate. In this section, we study the uniform pstability for the augmented system (4) with respect to initial data. Let Z := (Θ, Ω) andZ := (Θ,Ω) be two solutions to (4) corresponding to initial data (Θ 0 , Ω 0 ) and (Θ 0 ,Ω 0 ), respectively. For the uniform stability estimate, we introduce a metric which is equivalent to p -distance: for p ∈ [1, ∞) and two solutions Z = (Θ, Ω) andZ = (Θ,Ω), we define the distance as Next, we present a uniform p -stability of system (4) with respect to initial data as follows.
Definition 4.1. The system (4) is uniformly p -stable with respect to initial data if the following relation holds: For two solutions Z andZ to (4) with initial data Z 0 andZ 0 , respectively, there exists a positive constant G independent of t such that In the following lemma, we will derive differential inequalities for two subfunctionals Θ(t) −Θ(t) p and Ω(t) −Ω(t) p for p ∈ [1, ∞).
Proof. • Case A (Derivation of the first inequality (30) 1 ). Note that θ i −θ i satisfies

This yields
We multiply by p|θ i −θ i | p−1 on both sides, sum up the resulting relations with respect to i and apply Hölder's inequality to obtain This implies the desired estimate.
• Case B (Derivation of the first inequality (30) 2 ). Note that ω i −ω i satisfies where θ * ik is located between (θ k − θ i ) and (θ k −θ i ) by mean value theorem. By multiplying (ω i −ω i ) on both sides of (31), we have We use (32) and similar argument used in Proposition 1 to obtain By Hölder's inequality, we have Then, by using these estimation and relation (33), we obtain By applying the relation D(Ω) ≤ Ω p and Theorem 3.3, we attain the desired result.

SEUNG-YEAL HA, JEONGHO KIM, JINYEONG PARK AND XIONGTAO ZHANG
We combine Lemma 2.6 and Lemma 4.2 to obtain the uniform p -stability.
Theorem 4.3. Suppose that initial data and coupling strength satisfy the following relations: Then, for any two solutions (Θ, Ω) and (Θ,Ω), we have uniform p -stability estimate (29).
As a direct application of Theorem 4.3, we have the following corollary for the first-order Kuramoto model (1).
On the other hand, we also set initial frequencies: Then, we solve the second-order system (4) with initial data (34) and (35). It follows from the equivalence relation between KM (1) and AKM (4) in Theorem 2.5 and Theorem 4.3 that we have where Ω 0 = (ω 0 1 , · · · , ω 0 N ). We again use the relations (35) to find This yields Finally, we combine (36) and (37) to obtain the desired stability estimate.
5. Uniform mean-field limit from the AKM to kinetic equation. In this section, we present the uniform mean-field limit for the AKM in a measure theoretic framework. The limiting mean-field kinetic equation can be formally derived from the particle model (4) via the formal procedure of BBGKY hierarchy, and it can be rigorously justified using the standard empirical measure approximations and local-in-time stability estimates in Monge-Kantorovich distance which is equivalent to Wasserstein-1 distance in any finite time.
The formal BBGKY hierarchy procedure yields a formal mean-field limit of system (4) toward the mean-field kinetic equation as N → ∞. More precisely, let f = f (θ, ω, t) be the one-particle distribution function. Then, the kinetic equation reads as follows.
Recall that our main purpose of this section is to justify the rigorous transition from (4) to (38) in the mean-field limit (N → ∞).

5.1.
A measure theoretic framework. In this subsection, we briefly discuss some framework which embodies (4) and (38) in a common framework. For this, we first review concept of measure-valued solutions to (38).
Let P(T × R) be the set of all Radon probability measures with compact support on the phase space T × R, which can be understood as normalized nonnegative bounded linear functionals on C 0 (T × R). For a probability measure µ ∈ P(T × R), we use a standard duality relation: Next, we recall several definitions to be used later.
) be a measurevalued solution to (38) with initial data µ 0 ∈ P(T × R) if the following three assertions hold: 1. Total mass is normalized: µ t , 1 = 1. 2. µ is weakly continuous in t: Remark 4. Note that for a solution {(θ i , ω i )} to (4), the empirical measure is a measure-valued solution in the sense of Definition 5.1 to (38). Thus, ODE solution to (4) can be understood as a measure-valued solution for (38). Likewise, the classical solution for the kinetic AKM model (38) is also a measure-valued solution as well. Thus, we can treat the particle and kinetic AKM models in the same framework.
We now discuss how to measure the distance between the solutions of (4) and (38) by equipping a metric to the probability measure space P(T × R), and the concept of local-in-time mean-field limit. In fact, we can endow Wasserstein-p distance W p in the probability space P(T × R).
Definition 5.2. [33,40] 1. For p ∈ Z + , let P p (T × R) be a collection of all probability measures with finite p th moment: for some z 0 ∈ T × R µ, z − z 0 p p < +∞. Then, Wasserstein p-distance W p (µ, ν) is defined for any µ, ν ∈ P p (T × R) as where Γ(µ, ν) denotes the collection of all probability measures on T 2 × R 2 with marginals µ and ν. 2. If lim p→∞ W p exists, then we define W ∞ metric as the limit.
3. For any T ∈ (0, ∞], the kinetic equation (38) is derivable from the particle model (4) in [0, T ), or equivalent to say the mean-field limit from the particle system (4) to the kinetic equation (38), which is valid in [0, T ), if for every solution µ t of the kinetic equation (38) with initial data µ 0 , the following condition holds: for some p ∈ [1, ∞) and t ∈ [0, T ), where µ N t is a measure valued solution of the particle system (4) with initial data µ N 0 .
For later use, we quote two results on the approximation of a measure by empirical measures and mean-field limit in any finite time interval without proofs.

Proposition 2. [40]
For any given p ∈ [1, ∞) and µ ∈ P p (T × R) with compact support, there exists a sequence of empirical measures µ N ∈ P p (T × R) such that µ N has a common compact support with µ and lim N →+∞ Remark 5. The construction of the approximation can be followed by the method of Theorem 6.18 in the book [40] by finding a sequence of atomic measures N j=1 a j δ j with rational numbers a j such that N j=1 a j = 1.

5.2.
A uniform mean-field limit. In this subsection, we present a uniform meanfield limit to the kinetic equation (38). We basically follow the approach given in [21], Corollary 1 and Lemma 3.2.
On the other hand, since µ N t has a common compact support, we can view V p p as a test function. Then, due to Theorem 5.3, we have, This implies the weak convergence of µ N t to µ t . Thus, we can pass to the limit N → ∞ to (44) to obtain R 2d

5.3.
Complete synchronization estimate. In this subsection, we present an alternative approach for the complete synchronization estimate for (38) introduced in previous subsection. As mentioned in abstract, the kinetic equation (38) for the AKM is more suitable for the Lyapunov functional approach, compared to the kinetic Kuramoto equation for the KM with distributed natural frequencies. For simplicity of presentation, we suppress t-dependence in f : Lemma 5.4. Let f be a classical solution of (38) whose support is compact. Then, we have Proof. It directly comes from multiplying by 1 and ω to (38) and integrating the resulting relation over the phase and frequency space, hence we omit the detailed calculation.
Next, we discuss the derivation of the complete (frequency) synchronization estimate. For this, we use the Lyapunov functional defined as follows.
where ω c is the mean frequency which is constant due to Lemma 5.4. If complete synchronization occurs, it is natural to expect that frequency will converge to ω c , i.e., the Lyapunov functional Λ(f ) converges to 0. To show this, we will use the standard Lyapunov functional estimate on Λ(f ).
Theorem 5.5. Let f be a classical solution of (38) whose support is compact and initial datum f 0 satisfying where D 0 Θ is the diameter of support of f 0 projected to θ-space. Then, if the coupling strength κ is large enough, the Lyapunov functional Λ[f ] decays exponentially: Proof. It follows from Lemma 3.1 and condition for support of initial data that we have where D Θ (t) is the diameter of support of f (·, ·, t). Now we use (45) and use the periodicity of f in θ-variable to obtain Next, we estimate the terms I i separately.
• (estimate of I 1 ). By interchanging ω and ω * , we obtain where we use the condition |θ * − θ| ≤ D ∞ when θ and θ * are contained in support of f (·, ·, t), which is guaranteed by condition on support of initial data.
• (estimate of I 2 ). From the anti-symmetry of integrand, it is easy to see From the estimation of I i , i = 1, 2, we derive following differential inequality By using Grönwall's lemma, we can obtain the desired exponential decay.
6. Applications to the kinetic Kuramoto Model. In this section, we study the uniform mean-field limit of the Kuramoto model and as a direct application of previous results, we also show the existence of phase-locked states for the kinetic Kuramoto equation via uniform mean-field limit by lifting particle results to the kinetic regime. For the local-in-time stability and mean-field limit of the kinetic Kuramoto equation, we refer to [28].
Let f = f (θ, ν, t) be a one-oscillator probability density function for the ensemble of Kuramoto oscillators. Then, the dynamics of f is governed by the kinetic Kuramoto equation: Note that the probability density function g = g(ν) for natural frequencies appears as a ν-marginal density function of f : 2π 0 f (θ, ν, t)dθ = g(ν).
Theorem 6.2. Suppose that initial probability measure µ 0 ∈ P(T×R) and coupling strength satisfy .
Then, the following assertions hold. For p ∈ [1, ∞), second-order lifting of the Kuramoto model. For the corresponding kinetic equation, the complete frequency synchronization can be obtained via the Lyapunov functional measuring the dispersion of the frequency variations. Our proposed second-order model has a formal similarity to the Cucker-Smale flocking model. As long as the phase diameter is confined in a quarter arc, the flocking estimate, uniform p -stability and mean-field limit for the extended Kuramoto model can be analyzed using the similar techniques done for the Cucker-Smale model. As aforementioned in Introduction, the reason that we focus on the mean-field case (all-to-all couplings) is that we are interested in the complete synchronization estimate for the kinetic Kuramoto equation at the kinetic level without lifting particle results to the kinetic level. Other than this, some of the estimates studied in this paper can be extended to a more general setting. For example, we considered synchronization dynamics of Kuramoto oscillators on the complete network, but some estimates such as the uniform stability and synchronization estimates at the particle level can be certainly extended to the locally coupled Kuramoto oscillators over the general networks, say symmetric and connected networks. Moreover, nonlinear stability and instability of the incoherent state can be studied for the proposed kinetic equation.
We leave these interesting issues as a future work.