A hydrodynamic model for synchronization phenomena

TL;DR

This paper introduces a new hydrodynamic model derived from the Kuramoto model with inertia to study synchronization, revealing finite-time blow-up phenomena and providing analytical and numerical insights.

Contribution

The paper develops a novel pressureless Euler-based hydrodynamic model for synchronization, establishing well-posedness, synchronization conditions, and blow-up analysis, which extends classical Kuramoto models.

Findings

Finite-time blow-up phenomena observed in the model.

Synchronization estimates for identical natural frequencies.

Numerical simulations showing phase transitions and hysteresis.

Abstract

We present a new hydrodynamic model for synchronization phenomena which is a type of pressureless Euler system with nonlocal interaction forces. This system can be formally derived from the Kuramoto model with inertia, which is a classical model of interacting phase oscillators widely used to investigate synchronization phenomena, through a kinetic description under the mono-kinetic closure assumption. For the proposed system, we first establish local-in-time existence and uniqueness of classical solutions. For the case of identical natural frequencies, we provide synchronization estimates under suitable assumptions on the initial configurations. We also analyze critical thresholds leading to finite-time blow-up or global-in-time existence of classical solutions. In particular, our proposed model exhibits the finite-time blow-up phenomenon, which is not observed in the classical Kuramoto models, even with a smooth distribution function for natural frequencies. Finally, we numerically investigate synchronization, finite-time blow-up, phase transitions, and hysteresis phenomena.

Equations

$$\begin{array}{rl}& {\partial }_t\rho +{\partial }_{\theta }(\rho u)=0,(\theta ,\Omega )\in \mathbb{T}\times \mathbb{R},t>0,\\ & {\partial }_t(\rho u)+{\partial }_{\theta }(\rho u^2)=\frac{1}{m}\left(-\rho u+\rho \Omega +K\rho {\int }_{\mathbb{T}\times \mathbb{R}}\sin ({\theta }_{\ast }-\theta )\rho ({\theta }_{\ast },{\Omega }_{\ast },t)g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast }\right),\end{array}$$

$$(\rho ,u)(\theta ,\Omega ,0)=:({\rho }_0(\theta ,\Omega ),u_0(\theta ,\Omega )),(\theta ,\Omega )\in \mathbb{T}\times \mathbb{R}.$$

$$m{\ddot{\theta }}_i(t)+{\dot{\theta }}_i(t)={\Omega }_i+\frac{K}{N}\sum _{j=1}^N\sin ({\theta }_j(t)-{\theta }_i(t)),i=1,\cdots ,N,t>0.$$

$${\partial }_{t}f+{\partial }_{\theta }(\omega f)+{\partial }_{\omega }(F[f]f)=0,(\theta ,\omega ,\Omega )\in \mathbb{T}\times \mathbb{R}\times \mathbb{R},t>0,$$

$$F[f](\theta ,\omega ,\Omega ,t)=\frac{1}{m}\left(-\omega +\Omega +K{\int }_{\mathbb{T}\times \mathbb{R}\times \mathbb{R}}\sin ({\theta }_{\ast }-\theta )f({\theta }_{\ast },{\omega }_{\ast },{\Omega }_{\ast },t)g({\Omega }_{\ast })d{\theta }_{\ast }d{\omega }_{\ast }d{\Omega }_{\ast }\right).$$

$$\rho (\theta ,\Omega ,t)={\int }_{\mathbb{R}}fd\omega ,(\rho u)(\theta ,\Omega ,t)={\int }_{\mathbb{R}}\omega fd\omega ,$$

$$(\rho E)(\theta ,\Omega ,t)=\frac{1}{2}\rho u^2+\rho e,\text{mbox}where\rho e=\frac{1}{2}{\int }_{\mathbb{R}}|\omega -u{|}^{2}fd\omega .$$

$${\partial }_t\rho$$

$${\partial }_t(\rho u)$$

$$=-{\partial }_{\theta }(\rho u^2+p)-\frac{1}{m}\left(\rho u-\rho \Omega -\rho K{\int }_{\mathbb{T}\times \mathbb{R}}\sin ({\theta }_{\ast }-\theta )\rho ({\theta }_{\ast },{\Omega }_{\ast },t)g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast }\right),$$

$$p={\int }_{\mathbb{R}}|\omega -u{|}^{2}fd\omega ,$$

$${\int }_{\mathbb{R}}{\omega }^{2}fd\omega ={\int }_{\mathbb{R}}{u}^{2}fd\omega +{\int }_{\mathbb{R}}|\omega -u{|}^{2}fd\omega =\rho u^2+p.$$

$$\frac{1}{2}{\omega }^3=\frac{1}{2}(\omega -u+u{)}^3=\frac{1}{2}(\omega -u{)}^3+\frac{3}{2}(\omega -u{)}^{2}u+\frac{3}{2}(\omega -u)u^2+\frac{1}{2}{u}^3,$$

$${\partial }_t(\rho E)$$

$$=-{\partial }_{\theta }\left(q+\frac{3pu}{2}+\frac{1}{2}\rho u^3\right)-\frac{2}{m}\rho E+\frac{1}{m}\rho u\Omega$$

$$+\frac{K}{m}\rho u{\int }_{\mathbb{T}\times \mathbb{R}}\sin ({\theta }_{\ast }-\theta )\rho ({\theta }_{\ast },{\Omega }_{\ast },t)g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast }.$$

$$q=\frac{1}{2}{\int }_{\mathbb{R}}(\omega -u{)}^{3}fd\omega .$$

$$\begin{array}{rl}& {\partial }_t\rho +{\partial }_{\theta }(\rho u)=0,\\ & {\partial }_t(\rho u)+{\partial }_{\theta }(\rho u^2+p)\\ & =-\frac{\rho }{m}\left(u+\Omega -K{\int }_{\mathbb{T}\times \mathbb{R}}\sin ({\theta }_{\ast }-\theta )\rho ({\theta }_{\ast },{\Omega }_{\ast },t)g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast }\right),\\ & {\partial }_t(\rho E)+{\partial }_{\theta }\left(q+\frac{3pu}{2}+\frac{1}{2}\rho u^3\right)\\ & =-\frac{\rho }{m}\left(2E+u\Omega -Ku{\int }_{\mathbb{T}\times \mathbb{R}}\sin ({\theta }_{\ast }-\theta )\rho ({\theta }_{\ast },{\Omega }_{\ast },t)g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast }\right).\end{array}$$

$$f(\theta ,\omega ,\Omega ,t)\simeq \rho (\theta ,\Omega ,t){\delta }_{u(\theta ,\Omega ,t)}(\omega ).$$

$$f(\theta ,\omega ,\Omega ,t)\simeq \rho (\theta ,\Omega ,t)\exp \left(-\frac{(\omega -u(\theta ,\Omega ,t){)}^2}{2}\right),$$

$${\partial }_t\eta (\theta ,\Omega ,t)=u(\eta (\theta ,\Omega ,t),\Omega ,t)\text{mbox}with\eta (\theta ,\Omega ,0)=\theta .$$

$$h(\theta ,\Omega ,t):=\rho (\eta (\theta ,\Omega ,t),\Omega ,t)\text{mbox}andv(\theta ,\Omega ,t):=u(\eta (\theta ,\Omega ,t),\Omega ,t),$$

$$\begin{array}{rl}& h(\theta ,\Omega ,t){\partial }_{\theta }\eta (\theta ,\Omega ,t)={\rho }_0(\theta ,\Omega ),(\theta ,\Omega )\in {\mathcal{S}}_0,t>0,\\ & m{\partial }_{t}v(\theta ,\Omega ,t)+v(\theta ,\Omega ,t)\\ & =\Omega +K{\int }_{\mathbb{T}\times \mathbb{R}}\sin (\eta ({\theta }_{\ast },{\Omega }_{\ast },t)-\eta (\theta ,\Omega ,t)){\rho }_0({\theta }_{\ast },{\Omega }_{\ast })g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast },\end{array}$$

$$(h,v)(\theta ,\Omega ,0)=({\rho }_0,u_0)(\theta ,\Omega ),(\theta ,\Omega )\in {\mathcal{S}}_0.$$

$$\begin{array}{rl}& m{\partial }_{tt}\eta (\theta ,\Omega ,t)+{\partial }_t\eta (\theta ,\Omega ,t)\\ & =\Omega +K{\int }_{\mathbb{T}\times \mathbb{R}}\sin (\eta ({\theta }_{\ast },{\Omega }_{\ast },t)-\eta (\theta ,\Omega ,t)){\rho }_0({\theta }_{\ast },{\Omega }_{\ast })g({\Omega }_{\ast })d{\theta }_{\ast }d{\Omega }_{\ast },\end{array}$$

$$(\eta ,{\partial }_t\eta )(\theta ,\Omega ,0)=(\theta ,u_0(\theta ,\Omega )),(\theta ,\Omega )\in {\mathcal{S}}_0.$$

$$\begin{array}{rl}& {\partial }_t\rho +{\partial }_{\theta }(\rho u)=0,\theta \in \mathbb{T},t>0,\\ & {\partial }_t(\rho u)+{\partial }_{\theta }(\rho u^2)=\frac{1}{m}\left(-\rho u+\rho \Omega +K\rho {\int }_{\mathbb{T}}\sin ({\theta }_{\ast }-\theta )\rho ({\theta }_{\ast },t)d{\theta }_{\ast }\right).\end{array}$$

$$r(t)e^{\mathrm{i}\phi (t)}={\int }_{\mathbb{T}\times \mathbb{R}}{e}^{\mathrm{i}\theta }\rho (\theta ,\Omega ,t)g(\Omega )d\theta d\Omega ,$$

$$r(t)={\int }_{\mathbb{T}\times \mathbb{R}}\cos (\theta -\phi (t))\rho (\theta ,\Omega ,t)g(\Omega )d\theta d\Omega .$$

$${\partial }_t(\rho u)+{\partial }_{\theta }(\rho u^2)=-\frac{\rho }{m}\left(u-\Omega -Kr(t)\sin (\phi (t)-\theta )\right).$$

$$r(t)e^{\mathrm{i}\phi (t)}={\int }_{\mathbb{T}\times \mathbb{R}}{e}^{\mathrm{i}\eta (\theta ,\Omega ,t)}{\rho }_0(\theta ,\Omega )g(\Omega )d\theta d\Omega ,$$

$$m{\partial }_{t}v(\theta ,\Omega ,t)+v(\theta ,\Omega ,t)=\Omega +Kr(t)\sin (\phi (t)-\eta (\theta ,\Omega ,t)).$$

Full text

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 A priori estimates

4. 2.2 Auxiliary lemmas

5. 3 Local-in-time existence and uniqueness of classical solutions

6. 4 Synchronization estimates for identical oscillators

7. 4.1 Kinetic energy estimate

8. 4.2 Phase & velocity diameter estimates

9. 5 Critical thresholds phenomena

10. 6 Numerical experiments

11. 6.1 Numerical scheme

12. 6.2 Time evolutions of density and velocity

13. 6.3 Phase transitions & hysteresis phenomena

14. 7 Conclusion

15. A Proof of Theorem 5.1

1. Introduction

The phenomenon of collective synchronization exhibited by various biological systems is ubiquitous in nature, and it has been extensively studied in many different scientific disciplines such as applied mathematics, physics, biology, sociology, and control theory due to their biological and engineering applications, [1, 5, 17, 22, 27, 32, 34, 40, 44]. The mathematical treatment of synchronization phenomena was pioneered by Winfree [44] and Kuramoto [27]. Winfree first introduced a first-order model for collective synchronization of weakly coupled nonlinear oscillators. Subsequently, Kuramoto proposed a mathematically tractable model consisting of a population of coupled phase oscillators having natural frequencies extracted from a given distribution, and all of them are coupled by a mean-field interaction, sinusoidal coupling. The Kuramoto model contains all the main features of interest. In particular, the Kuramoto model displays a phase transition between coherent and incoherent states: the oscillators rotate on a circle incoherently when the coupling strength is weak enough, while the collective synchronization occurs when the coupling strength is beyond some threshold. Since then, the Kuramoto model has become a paradigmatic model for synchronization phenomena.

There already exist various extensions, such as additive/multiplicative noises, time-delayed coupling, inertia, frustration, and networks, extensively explored in [4, 7, 11, 12, 14, 15, 22, 23, 24, 37, 38]. However, to the best of our knowledge, there is no available literature on hydrodynamic models for synchronization phenomena. In the current work, we present a new hydrodynamic model, which is a pressureless Euler-type system, for the synchronization phenomena. More specifically, let $\rho=\rho(\theta,\Omega,t)$ and $u=u(\theta,\Omega,t)$ be the density and velocity functions of Kuramoto oscillators, respectively, in $\theta\in\mathbb{T}:=\mathbb{R}/(2\pi\mathbb{Z})$ with a natural frequency $\Omega$ extracted from a given distribution function $g=g(\Omega)$ at time $t>0$. Then our main system is given by

[TABLE]

with the initial data

[TABLE]

Here $m>0$ and $K>0$ denote the strength of inertia and coupling strength, respectively.

The system (1.1) can be formally derived from the second-order system of ordinary differential equations for synchronization, called the Kuramoto model with inertia, through a kinetic description under a mono-kinetic closure assumption. More precisely, our starting point is the $N$-particle Kuramoto oscillators with inertia. Let $\theta_{i}\in\mathbb{R}$ be the phase of the $i$-th oscillator with the natural frequency $\Omega_{i}$. Then the dynamics of second-order Kuramoto oscillators is governed by the following system:

[TABLE]

The particle model (1.3) is introduced in [22] as a phenomenological model to describe the slow relaxation in the synchronization process in certain biological systems, e.g., fireflies of the Pteroptyx malaccae. Note that the classical Kuramoto model can be simply obtained by disregarding the inertial effect, i.e., setting $m=0$. A different set of applications of the second-order phase model (1.3) includes power grids, superconducting Josephson junction arrays, and explosive synchronization [15, 19, 20, 25, 41, 42, 43]. Furthermore, it is known that the model (1.3) exhibits rich phenomena such as the discontinuous phase transition and hysteretic dynamics [4, 37, 38]. For mathematical results on (1.3), we refer to [11, 13, 14, 16, 20].

On the other hand, when the number of oscillators $N$ is very large, the microscopic description (1.3) is computationally complicated, and thus understanding how this complexity can be reduced is an important issue. The classical strategy to reduce this complexity is to derive a mesoscopic description, i.e., continuum model of the dynamics, by introducing a distribution function. Let $f=f(\theta,\omega,\Omega,t)$ be the one-oscillator distribution function on the space $\mathbb{T}\times\mathbb{R}$ with the natural frequency $\Omega$ at time $t$ and satisfy the normalized condition $\int_{\mathbb{T}\times\mathbb{R}}f(\theta,\omega,\Omega,t)\,d\theta d\omega=1$. At the formal level, we can expect that as the number of oscillators $N$ goes to infinity, the $N$-particle system (1.3) will be replaced by the following Vlasov-type equation:

[TABLE]

where the interaction term $F[f]$ is given by

[TABLE]

The kinetic equation (1.4) is often used in the physics literature to study the phase transition phenomena [2, 3]. The rigorous derivation of the equation (1.4) from (1.3) is established in [12]. The global-in-time existence of measure-valued solutions and its long-time behavior are also studied in [12]. We refer to [18, 31, 33] for the rigorous derivation of kinetic equations.

Note that the mesoscopic description model (1.4) is posed in $3+1$ dimensions, thus obtaining a numerical solution of (1.4) is computationally expensive. For this reason, we next derive a macroscopic description model from (1.4) by taking moments together with zero temperature closure or mono-kinetic assumption for the local hydrodynamic solutions. In this way, we can remove $\omega$-variable in solutions. For this, we first set local density $\rho$, moment $\rho u$, and energy $\rho E$, which is the sum of kinetic and internal energies:

[TABLE]

and

[TABLE]

Then straightforward computations yield

[TABLE]

where we denote by $p$ the pressure given by

[TABLE]

and we used

[TABLE]

Since

[TABLE]

we also find

[TABLE]

Here we denote by $q$ the heat flux given by

[TABLE]

By collecting the above observations, we have the following local conservation laws:

[TABLE]

In order to close the local conservation laws (1.5), we use the mono-kinetic closure assumption:

[TABLE]

Then this reduces to our main system (1.1). Although the mono-kinetic assumption is not fully justified, it is known that the hydrodynamic system derived gives quantitative results comparable to the particle simulations, see [8, 9, 10].

Remark 1.1**.**

We can also employ another closure assumption, for instance, a local Maxwellian-type ansatz,

[TABLE]

then we derive the isothermal Euler-type equations from (1.5).

In the current work, we first establish local-in-time existence and uniqueness of classical solutions to the system (1.1). For this, we consider the moving domain problem and reformulate our main system (1.1) into the Lagrangian coordinate. To be more precise, let us define the characteristic flow $\eta(\theta,\Omega,t)$ by

[TABLE]

Set

[TABLE]

and let us denote by the time-varying set $\mathcal{S}_{t}:=\{(\theta,\Omega)\in\mathbb{T}\times\mathbb{R}:\rho(\theta,\Omega,t)\neq 0\}$ for given initially bounded open set $\mathcal{S}_{0}$. Using these newly defined notations, we can rewrite the system (1.1) along the characteristic flow given in (1.6) as

[TABLE]

with the initial data

[TABLE]

For the system (1.7), we introduce a weighted Sobolev space $H^{s}_{g}$ by the distribution function $g$ and construct a unique $H^{s}_{g}$-solution. This newly defined solution space together with our careful analysis allows us to apply directly our strategy for the identical case, i.e., the distribution function $g$ is the form of the Dirac measure on $\mathbb{R}$ giving unit mass to the point $\Omega_{0}$, $g(\Omega)=\delta_{\Omega_{0}}(\Omega)$ for some $\Omega_{0}\in\mathbb{R}$. We construct the approximated solutions and provide that they are Cauchy sequences in the proposed weighted Sobolev spaces by obtaining uniform bound estimates of approximated solutions. We then show that the limiting functions are solutions to (1.7). The details of proof are discussed in Section 3. It is worth noticing that our system is a type of pressureless Euler equations with nonlocal forces, and it is well known that the pressureless Euler system may develop a singularity such as a $\delta$-shock in finite-time, i.e., fail to admit a global classical solution, no matter how smooth the initial data are. This is one of the main difficulties in analyzing the Euler equations. However, for the identical case, in general the case where the distribution function $g$ is a sum of Dirac measures, we expect that the density $\rho$ converges toward a Dirac measure, see Remark 4.4. This infers that the existence time of $H^{s}_{g}$-solutions cannot be infinity for general $g$.

Remark 1.2**.**

Since $\partial_{t}v=\partial_{tt}\eta$, we can also rewrite the momentum equation in (1.7) as

[TABLE]

subject to the initial data

[TABLE]

Note that the equation (1.9) is a closed equation, i.e., the continuity equation in (1.7) is decoupled from the equation for $v$.

After we construct the local-in-time existence of classical solutions, we discuss the synchronization estimate for the case of identical oscillators in Section 4. In this case, upon rotating frame if necessary, we may assume that $g(\Omega)=\delta_{0}(\Omega)$ and the system (1.1) reduces to

[TABLE]

For the system (1.10), we present two different methods for the synchronization estimates. Inspired by a recent work [13], in Section 4.1, we propose a strategy based on kinetic energy combined with order parameter $r$ estimates, where $r$ is defined by

[TABLE]

i.e.,

[TABLE]

Here $\varphi$ represents the average phase associated to the system (1.1). Note that the order parameter $r$ is employed to measure the phase transition from incoherent to coherent states mentioned above. To be more specific, $r^{\infty}:=\lim_{t\to\infty}r(t)$ as a function of the coupling strength $K$, i.e., $r^{\infty}=r^{\infty}(K)$, changes from zero (incoherent state or disordered state) to a non-zero value (coherent state or ordered state) when the coupling strength $K$ exceeds a critical value $K_{c}$. It is known that the critical coupling strength is $K_{c}=2/(\pi g(0))$ for the classical Kuramoto model (1.3) with $m=0$ when $g$ is unimodal and symmetric about $\Omega=0$. Our synchronization estimate provides the convergences of the velocity toward the mean velocity and the order parameter to some positive constant under general assumptions on the initial configurations, this strategy can be applied for the case $\mathcal{S}_{0}=\mathbb{T}\times\mathbb{R}$, see Remark 3.3 and Theorem 4.1. However, it does not give any information about the decay rate of convergence and the limit profiles $\rho^{\infty}:=\lim_{t\to+\infty}\rho$.

Remark 1.3**.**

Using newly defined notations in (1.11), we can rewrite the momentum equation in (1.1) as

[TABLE]

For the system (1.7), the order parameter $r$ and the average phase $\varphi$ can be expressed by

[TABLE]

and the equation of velocity in (1.7) can be rewritten as

[TABLE]

In order to complement the drawbacks of the previous strategy, in Section 4.2, we provide a second-order Grönwall-type inequality estimates on phase and velocity diameters for the synchronization. Note that the equation for velocity in (1.7) is an integro-differential equation, and the equation (1.9) resembles the particle Kuramoto model with inertia (1.3). In view of this fact, we use the idea of [11] and estimate the phase and velocity diameters to show the exponential synchronization behavior under certain assumptions on the initial data. Although this approach requires more restricted class of initial data than that in the previous approach, it gives decay rates of convergences and shows the limit profiles $\rho^{\infty}$ is the form of Dirac measure, see Remark 4.4.

In Section 5, we show that our main system (1.1) exhibits critical threshold phenomena. For this, we first study the local-in-time well-posedness of the system (1.1) if the initial data are sufficiently regular and the initial density has no vacuum. Then, we analyze critical thresholds determining regions of initial conditions for global-in-time existence and finite-time blow-up of solutions to the system (1.1). More precisely, we provide thresholds between the supercritical region with finite-time breakdown and the subcritical region with global-in-time existence of the classical solutions. The critical threshold phenomenon for Eulerian dynamics is studied in [21, 30, 36] for Euler-Poisson equations and [6, 35] for pressureless Euler equations with nonlocal velocity alignment forces. We want to emphasize that the finite-time blow-up of solutions cannot be observed in the Kuramoto model with inertia at both microscopic level (1.3) and mesoscopic level (1.4). As mentioned above, it is an important issue for the global existence of solutions to the Euler-type equations how to prevent the formation of singularity. However, this implies that our hydrodynamic model (1.1) may describe the finite-time synchronization phenomena [29, 45] commonly found in some natural networks, see also Remark 5.4. We investigate the supercritical region for the system (1.1) so that the classical solution will blow up in finite time if its initial data belong to that region. On the other hand, we show that if the initial data is in the subcritical region, then the initial regularity of solutions is preserved in time. These results are stated in Theorem 5.2.

Numerical experiments validating our theoretical results and giving further insights are presented in Section 6. We employ a finite-volume type scheme for the numerical simulations. We use initial data and parameters lying in sub/supercritical regions to illustrate the time evolution of solutions $\rho$ and $u$. The numerical simulations show that the finite-time blow-up of solutions may imply that the formation of finite-time synchronization, see Figures 1, 2, 4, and 5. It is also very interesting that our main system (1.1) also exhibits the hysteresis phenomenon. Depending on the strength of $m$, we show different types of phase transitions of the order parameter $r^{\infty}(K)$, see Figure 6. We would like to emphasize that our hydrodynamic model (1.1) is much more efficient than the $N$-particle system (1.3) in terms of computational cost when $N$ is large. We finally summarize our main results and report future research directions in the last section.

Before closing this section, we introduce several notations used throughout the paper. For a function $f(\theta,\Omega)$, $\|f\|_{L^{p}}$ denotes the usual $L^{p}(\mathbb{T}\times\mathbb{R})$-norm. $L^{p}_{g}(\mathbb{T}\times\mathbb{R})$ represents the space of weighted measurable functions whose $p$-th powers weighted by $g=g(\Omega)$ are integrable on $\mathbb{T}\times\mathbb{R}$, with the norm

[TABLE]

$f\lesssim g$ implies that there exists a positive constant $C>0$ such that $f\leq Cg$. We also denote by $C$ a generic, not necessarily identical, positive constant. For any nonnegative integer $s$, $H^{s}$ and $H^{s}_{g}$ denote the $s$-th order $L^{2}$ and $L^{2}_{g}$ Sobolev spaces, respectively. $\mathcal{C}^{s}([0,T];E)$ is the set of $s$ times continuously differentiable functions from an interval $[0,T]\subset\mathbb{R}$ into a Banach space $E$, and $L^{p}(0,T;E)$ is the set of $L^{p}$ functions from an interval $(0,T)$ to a Banach space $E$.

Throughout this paper, we assume that the distribution function $g=g(\Omega)$ satisfies

[TABLE]

Note that the case of identical oscillators, $g(\Omega)=\delta_{\Omega_{0}}$ for some $\Omega_{0}\in\mathbb{R}$, satisfies the above conditions (1.13). Furthermore, we assume that the initial density $\rho_{0}$ satisfies

[TABLE]

2. Preliminaries

In this section, we present a priori energy estimates and some useful lemmas which will be frequently used later.

2.1. A priori estimates

We first provide a priori energy estimates for the system (1.1).

Lemma 2.1**.**

Let $(\rho,u)$ be a global classical solution to the system (1.1). Then we have

[TABLE]

Proof.

(i) It clearly follows from the continuity equation in (1.1).

(ii) Multiplying the momentum equation in (1.1) by $g(\Omega)$ and integrating the resulting relation with respect to $\theta$ and $\Omega$, we find

[TABLE]

since $\sin(-\theta)=-\sin(\theta)$ for $\theta\in\mathbb{T}$. On the other hand, it again follows from the continuity equation in (1.1) that

[TABLE]

This together with (2.1) asserts (ii).

(iii) A straightforward computation yields

[TABLE]

Then, we use the momentum equation in (1.1) to find

[TABLE]

By exchanging $(\theta,\Omega)\leftrightarrow(\theta_{*},\Omega_{*})$, we can rewrite the last term on the right hand side of the above equation as

[TABLE]

Combining (2.2) and (2.3), we have the desired result. ∎

Remark 2.1**.**

From the continuity equation in (1.1), we easily get

[TABLE]

due to (1.14). Then this together with Lemma 2.1 (i) and (1.13) yields

[TABLE]

Remark 2.2**.**

It follows from Lemma 2.1 (ii) that

[TABLE]

In particular, if we consider the case of identical oscillators, i.e., $g(\Omega)=\delta_{0}$, upon shifting if necessary, then the above estimate gives the exponential decay of the momentum:

[TABLE]

2.2. Auxiliary lemmas

In this part, we present several useful lemmas that will be used out later. We first provide the exponential decay estimates for the nonnegative functions satisfying the following second-order differential inequality:

[TABLE]

where $a>0,b$ and $c$ are constants. We recall [11, Lemma 3.1] the following inequalities.

Lemma 2.2**.**

Let $x=x(t)$ be a nonnegative $\mathcal{C}^{2}$-function satisfying the differential inequality (2.4).

• (i)

*If $b^{2}-4ac>0$, then we have *

[TABLE]

where decay exponents $\alpha_{1}$ and $\alpha_{2}$ are given by

[TABLE]

respectively.

• (ii)

*If $b^{2}-4ac\leq 0$, then we have *

[TABLE]

We next provide the following simple lemma without the proof.

Lemma 2.3**.**

Suppose that a real-valued function $f:[0,\infty)\to\mathbb{R}$ is uniformly continuous and satisfies

[TABLE]

Then, $f$ tends to zero as $t\to\infty$:

[TABLE]

We also present a decay estimate for some differential equation, the proof of which can be found in [13, Lemma 4.1].

Lemma 2.4**.**

Let $y=y(t)$ be a nonnegative $\mathcal{C}^{1}$-function satisfying

[TABLE]

where $\alpha>0$ and $\beta$ is a bounded continuous function decaying to zero as t goes to infinity. Then $y$ satisfies

[TABLE]

In particular, $y$ tends to zero as $t$ goes to infinity.

We finally recall from [26, Lemma 2.4] the following Sobolev inequality.

Lemma 2.5**.**

Let $k\geq 1$. For $f\in(H^{k}\cap L^{\infty})(\mathbb{T})$, let $p\in[1,\infty]$, and $h\in\mathcal{C}^{k}(B(0,\|f\|_{L^{\infty}}))$ where $B(0,R)$ denotes the ball of radius $R>0$ centered at the origin in $\mathbb{R}$. Then there exists a positive constant $C=C(k,p,h)$ such that

[TABLE]

3. Local-in-time existence and uniqueness of classical solutions

In this section, we present the local-in-time well-posedness of the Lagrangian system (1.7). More precisely, we show the local-in-time existence and uniqueness of $H^{s}_{g}$-solutions with $s\geq 1$ to the system (1.7).

Theorem 3.1**.**

Let $s\geq 1$. Suppose that $(\rho_{0},u_{0})\in H_{g}^{s}(\mathcal{S}_{0})\times H_{g}^{s+1}(\mathcal{S}_{0})$. Then for any constants $0<M<M^{\prime}$ there exists $T_{0}>0$ depending only on $M$ and $M^{\prime}$, such that if $\|u_{0}\|_{H_{g}^{s+1}}<M$, then the system (1.7)-(1.8) has a unique solution

[TABLE]

satisfying

[TABLE]

Remark 3.1**.**

Almost the same argument as above can be applied to the case of identical oscillators, i.e., $g(\Omega)=\delta_{\Omega_{0}}(\Omega)$. In this case, the weighted spaces $L^{p}_{g}(\mathbb{T}\times\mathbb{R})$ and $H^{s}_{g}(\mathbb{T}\times\mathbb{R})$ reduce to the spaces $L^{p}(\mathbb{T})$ and $H^{s}(\mathbb{T})$, respectively.

Proof of Theorem 3.1.

Although the proof is similar to that of [9, Theorem A.1], for the completeness of our paper, we provide the details of it. We first approximate solutions of the system (1.7) by the sequence $(\eta^{n},v^{n})$ which is the solution of the following integro-differential system:

[TABLE]

with the initial data and the first iteration step defined by

[TABLE]

and

[TABLE]

Here the interaction term $f^{n+1}$ is given by

[TABLE]

From now on, for the notational simplicity, we suppress the $\theta$- and $\Omega$-dependences of the variables and domain if the context is clear.

$\bullet$ (Step 1: Uniform bounds): We claim that there exists $T_{0}>0$ such that

[TABLE]

(Proof of claim): We use an induction argument. In the first iteration step, we find that

[TABLE]

Let us assume that

[TABLE]

for some $T>0.$ Then, we check that the linear approximations $(\eta^{n+1},v^{n+1})$ from the system $\eqref{ap}$ are well-defined, and they satisfy

[TABLE]

We begin by estimating $\eta^{n+1}$. It follows from the equation of $\eta^{n+1}$ in (LABEL:ap) that

[TABLE]

for $k\geq 1$, where $\delta_{k,1}$ denotes Kronecker delta, i.e., $\delta_{k,1}=1$ for $k=1$ and $\delta_{k,1}=0$ otherwise. The equation of $v^{n+1}$ in (LABEL:ap) gives

[TABLE]

Since

[TABLE]

we find

[TABLE]

Thus we obtain

[TABLE]

for some $C>0$, due to (1.13). For $1\leq k\leq s+1$, by taking $\partial_{\theta}^{k}$ to (3.2), we get

[TABLE]

We next use the facts $\sin\in\mathcal{C}^{\infty}(\mathbb{T})$ and $\|\sin\|_{L^{\infty}}\leq 1$ with Lemma 2.5 to estimate

[TABLE]

for some $C>0$. This yields

[TABLE]

We then use

[TABLE]

to have

[TABLE]

for some $C>0$. Combing the estimate above with (3.3) asserts

[TABLE]

Note that at $t=0$, the right hand side of (3.4) is $\|u_{0}\|_{H^{s+1}}\leq M<M^{\prime}$, thus there exists a small time $0<T_{0}<T$ such that

[TABLE]

$\bullet$ (Step 2: Cauchy estimates): For notational simplicity, we set

[TABLE]

Then we easily find

[TABLE]

thus

[TABLE]

and also we have

[TABLE]

Note that

[TABLE]

where we used Hölder’s inequality and $\|\rho_{0}\|_{L^{1}_{g}}=1$. Thus we get

[TABLE]

Then (3.6) and (3.7) yield

[TABLE]

By setting $\Delta^{n+1}_{\eta,v}(t):=\|\eta^{n+1,n}(\cdot,\cdot,t)\|_{L^{2}_{g}}+\|v^{n+1,n}(\cdot,\cdot,t)\|_{L^{2}_{g}}$, and then we combine (3.5) and (3.8) to get

[TABLE]

which implies that $\Delta^{n+1}_{\eta,v}(t)\lesssim T_{0}^{n+1}/(n+1)!$. Thus, we conclude that $(\eta^{n}(\theta,\Omega,t),v^{n}(\theta,\Omega,t))$ is a Cauchy sequence in $\mathcal{C}([0,T_{0}];L^{2}_{g}(\mathcal{S}_{0}))\times\mathcal{C}([0,T_{0}];L^{2}_{g}(\mathcal{S}_{0}))$.

$\bullet$ (Step 3: Regularity of limit functions): It follows from Step 2 that there exist limit functions $\eta$ and $v$ such that

[TABLE]

Interpolating this with the uniform bound estimates in Step 1, we obtain

[TABLE]

Now, we claim that

[TABLE]

Note that $v\in\mathcal{C}([0,T_{0}];H^{s+1}_{g}(\mathcal{S}_{0}))$ implies $\eta\in\mathcal{C}^{1}([0,T_{0}];H^{s+1}_{g}(\mathcal{S}_{0}))$ thanks to (3.9) and (LABEL:ap), so it suffices to show that $v\in\mathcal{C}([0,T_{0}];H^{s+1}_{g}(\mathcal{S}_{0}))$. It follows from Step 1 that there exists a subsequence $(v^{n_{k}})$ such that

[TABLE]

for some $\tilde{v}\in L^{\infty}([0,T_{0}];H^{s+1}_{g}(\mathcal{S}_{0})).$ On the other hand, we also have

[TABLE]

Thus, we conclude that

[TABLE]

which yields

[TABLE]

Next, we show that

[TABLE]

for $t_{0}\in[0,T_{0}]$. Without loss of generality, we may assume that $t_{0}=0$. Since $\|\cdot\|_{H^{s+1}}$ is weakly lower semicontinuous,

[TABLE]

To show the weak continuity, let $\{t_{k}\}$ be a sequence such that $t_{k}\to 0$, $\{t_{k}\}\subset[0,T]$. For this sequence, we already have $\lim_{k\to\infty}\|v(t_{k})-u_{0}\|_{H^{s}_{g}}=0$ and $\|u_{0}\|_{H^{s+1}_{g}}\leq M^{\prime}$. On the other hand, it is easy to see that (3.4) yields

[TABLE]

From (3.11) and (3.12), we get

[TABLE]

This together with the weak convergence (3.10), we conclude that

[TABLE]

By considering the time-reversed problem, i.e., $t\mapsto-t$, we can also obtain the left continuity in the same way.

$\bullet$ (Step 4: Existence): In Step 3, we obtained (3.9), and this implies that the limit functions $\eta$ and $v$ are the solutions of (1.6)-(1.7) in the sense of distributions. Moreover, we also have $h\in\mathcal{C}([0,T_{0}];H^{s}_{g}(\mathcal{S}_{0}))$ since $\rho_{0}\in H^{s}_{g}(\mathcal{S}_{0})$ and $\eta\in\mathcal{C}([0,T_{0}];H^{s+1}_{g}(\mathcal{S}_{0}))$.

$\bullet$ (Step 5: Uniqueness): Let $(h,\eta,v)$ and $(\tilde{h},\tilde{\eta},\tilde{v})$ be the two classical solutions with the same initial data $(\rho_{0},u_{0})$. Let $\eta$ and $\tilde{\eta}$ be the trajectories with respect to $v$ and $\tilde{v}$, respectively, or, equivalently,

[TABLE]

Then, by the similar arguments as in Step 2, we obtain the Grönwall’s inequality:

[TABLE]

which yields

[TABLE]

Hence, one can easily check that

[TABLE]

using the similar argument in Step 3. Finally, this yields

[TABLE]

∎

Remark 3.2**.**

Theorem 3.1 gives the local-in-time regularity of solutions for the Cauchy problem in the Lagrangian coordinates (1.7). In order to go back to the Cauchy problem in the Eulerian coordinates (1.1), we need to show that the characteristic flow $(\theta,\Omega)\mapsto(\eta(\theta,\Omega,t),\Omega)$ defined in (1.6) is a diffeomorphism for all $t\in(0,T_{0})$ for some $T_{0}>0$. For this, it suffices to show that $\det(\nabla_{(\theta,\Omega)}(\eta(\theta,\Omega,t),\Omega))=\partial_{\theta}\eta(\theta,\Omega,t)>0$ for all $(\theta,\Omega,t)\in\mathbb{T}\times\mathbb{R}\times(0,T_{0})$. However, it is unclear how to show this. On the other hand, for the identical oscillators, i.e., the system (1.10), it follows from (1.6) that

[TABLE]

Then, by Theorem 3.1, we have

[TABLE]

Thus, by choosing small enough $T_{0}>0$, we obtain

[TABLE]

This together with Theorem 3.1 concludes the local-in-time existence and uniqueness of solutions $(\rho,u)$ to the system (1.10) such that $(\rho,u)\in\mathcal{C}([0,T_{0}];H^{s}(\mathcal{S}_{t}))\times\mathcal{C}([0,T_{0}];H^{s+1}(\mathcal{S}_{t}))$ for some $T_{0}>0$ with $(\rho_{0},u_{0})\in H^{s}(\mathcal{S}_{0})\times H^{s+1}(\mathcal{S}_{0})$. Here $\mathcal{S}_{t}=\{\theta\in\mathbb{T}:\rho(\theta,t)\neq 0\}$.

Remark 3.3**.**

We can also directly study the existence of classical solutions to the system (1.1) without introducing the Lagrangian formulation (1.7) under the assumption that the initial density $\rho>0$ in $\mathbb{T}\times\mathbb{R}$, see Theorem 5.1. In this case, however, we cannot use the strategy in Section 4.2 for the synchronization estimate for the case of identical oscillators since it requires that the diameter of support of the initial density in phase is less than $\pi$, see Lemma 4.3.

4. Synchronization estimates for identical oscillators

In this section, we provide synchronization estimates for identical oscillators, i.e., $g(\Omega)=\delta_{\Omega_{0}}(\Omega)$ for some $\Omega_{0}\in\mathbb{R}$. Without loss of generality, upon rotating frame if necessary, we may assume $\Omega_{0}=0$. Note that in this case the system (1.1) reduces to

[TABLE]

In the Lagrangian formulation, it is given by

[TABLE]

As mentioned in Introduction, we propose two different types of strategies for the synchronization estimates in the following two subsections.

4.1. Kinetic energy estimate

We introduce the mean velocity and mean phase:

[TABLE]

respectively, and the kinetic and potential energy functions $\mathcal{E}_{k}$ and $\mathcal{E}_{p}$:

[TABLE]

Note that the quantities above can be reformulated in the Eulerian coordinate as follows: the mean velocity and mean phase are given by

[TABLE]

and the corresponding kinetic and potential energy functions are

[TABLE]

respectively.

Lemma 4.1**.**

Let $(\eta,h,v)$ be a global solution to the system (4.1). Then we have the following estimates.

• (i)

Mean velocity estimate:

[TABLE]

• (ii)

Mean phase estimate:

[TABLE]

• (iii)

Total energy estimate:

[TABLE]

Proof.

(i) Note that

[TABLE]

Clearly, the second term vanishes and the desired estimate follows.

(ii) It directly follows from (i).

(iii) It is clear that

[TABLE]

On the other hand, we use the equation for $v$ in (4.1) to find

[TABLE]

where $I_{i},i=1,2$ can be estimated as follows.

[TABLE]

Combining the above estimates concludes the desired result. ∎

Remark 4.1**.**

The function $\mathcal{E}_{p}$ can be rewritten in terms of the order parameter defined in (1.11):

[TABLE]

We now state our main results on the decay of kinetic energy $\mathcal{E}_{k}$ and the convergence of order parameter $r$.

Theorem 4.1**.**

Let $(\eta,h,v)$ be a global solution to the system (4.1). Then we have

[TABLE]

Proof.

$\bullet$ (Decay of kinetic energy): In view of Lemmas 2.3 and 4.1, it suffices to show that the kinetic energy function $\mathcal{E}_{k}(t)$ is uniformly continuous since $\mathcal{E}_{k}\in L^{1}(0,\infty)$. Note that the kinetic energy $\mathcal{E}_{k}$ satisfies

[TABLE]

We then estimate $I_{2}$ as

[TABLE]

This yields

[TABLE]

Since Lemma 4.1 (iii) implies $\mathcal{E}_{k}(t)\leq\mathcal{E}_{k}(0)+\mathcal{E}_{p}(0)$, we obtain

[TABLE]

$\bullet$ (Convergence of order parameter): It follows from Lemma 4.1 that

[TABLE]

Since $t\mapsto\int_{0}^{t}\mathcal{E}_{k}(s)\,ds$ is increasing and bounded by $(m/2)\mathcal{E}_{k}(0)+(K/4)(1-r_{0}^{2})$ from above, it converges. On the other hand, $\mathcal{E}_{k}(t)$ decays to zero as $t\to\infty$, and thus we get

[TABLE]

due to $r\geq 0$. This completes the proof. ∎

In the rest of this subsection, we further study the time evolution of solutions to the system (4.1). For this, we set

[TABLE]

We then show the convergence of $\mathcal{L}(t)$ to zero as $t$ goes to infinity in the proposition below.

Proposition 4.1**.**

Let $(\eta,h,v)$ be a global solution to the system (4.1). Then we have the following assertions.

• (i)

If

[TABLE]

then we have

[TABLE]

• (ii)

The function $\mathcal{L}(t)$ decays to zero as $t\to\infty$:

[TABLE]

In particular, if (4.2) holds, we have

[TABLE]

Proof.

(i) It is clear from Theorem 4.1 that $\lim_{t\to\infty}r(t)>0$. Suppose that there exists $t_{0}\in(0,\infty)$ such that

[TABLE]

On the other hand, it follows from Lemma 4.1 (iii) and Remark 4.1 that

[TABLE]

and taking the limit $t\to t_{0}-$ gives the following contradiction.

[TABLE]

Thus we have

[TABLE]

(ii) First, we notice that the following identity holds:

[TABLE]

Taking the time derivative to $\mathcal{L}$, we find

[TABLE]

Using (1.12), we can easily estimate $J_{1}$ as

[TABLE]

For the estimate of $J_{2}$, we use (4.3) to find

[TABLE]

Note that

[TABLE]

Thus we obtain

[TABLE]

for some $\varepsilon>0$ which will be determined later.

Combining all of the above estimates, we have

[TABLE]

We now choose $\varepsilon>0$ small enough such that $1-\sqrt{2}K\varepsilon>0$ and use Lemma 2.4 and Thoerem 4.1 to obtain

[TABLE]

Finally, since

[TABLE]

it follows from (4.5), Thoerem 4.1, and Lemma 4.1 (i) that

[TABLE]

This together with the result (i) concludes our desired result. ∎

Remark 4.2**.**

Since $\mathcal{E}_{k}\in L^{1}(0,\infty)$, we find from (4.4) that $\mathcal{L}\in L^{1}(0,\infty)$. This further yields

[TABLE]

Moreover, if the initial data satisfy (4.2), we have

[TABLE]

4.2. Phase & velocity diameter estimates

In this part, we provide phase and velocity diameter estimates showing the exponential synchronization behavior under certain assumptions on the initial configurations. For this, we introduce the phase and velocity diameter functions as follows.

[TABLE]

For the synchronization estimates, we derive Grönwall-type differential inequalities for $d_{\eta}$ and $d_{v}$. We first show differentiability of these functions in the lemma below.

Lemma 4.2**.**

Let $T>0$ and $(\eta,h,v)$ be a solution to the system (4.1) on the interval $[0,T]$. Then, there exists a $T^{*}\in(0,T]$ such that

• (i)

$\partial_{\theta}\eta(\theta,t)>0~{}\mbox{for all}~{}(\theta,t)\in\mathbb{T}\times[0,T^{*}),$**

• (ii)

$d_{\eta}^{\prime}(t)~{}\mbox{and}~{}d_{\eta}^{\prime\prime}(t)~{}\mbox{exist on the interval}~{}[0,T^{*}).$**

Proof.

(i) Note that $\eta$ defined in (1.6) satisfies

[TABLE]

Thus, we get

[TABLE]

Thus, the set defined by

[TABLE]

is nonempty, and thus the assertion (i) is obtained for $T^{*}:=\sup\mathcal{T}$.

(ii) In view of (i), we find that there is no intersection between the characteristic curves starting from different points $\theta$ and $\theta_{*}$ on $[0,T^{*})$. Accordingly, indices $M(t)$ and $m(t)$ which give $d_{\eta}(t)=\eta(\theta_{M(t)},t)-\eta(\theta_{m(t)},t)$ stay fixed on that time interval, i.e., the indices $M(t)$ and $m(t)$ are constants on $[0,T^{*})$. This yields that $d_{\eta}$ and $d_{\eta}^{\prime}$ are differentiable. ∎

We now set $C_{0}$ and $D_{0}$ as follows.

[TABLE]

Lemma 4.3**.**

(Uniform boundedness of phase diameter) Let $T>0$ and $(\eta,h,v)$ be a solution to the system (4.1) on the interval $[0,T]$. Suppose that the initial data satisfy $0<C_{0}<\pi.$ Then, we have

[TABLE]

where $T^{*}$ appeared in Lemma 4.2 and $C_{0}>0$ is independent of $t$.

Proof.

Suppose that the phase diameter satisfies $d_{\eta}(t)<\pi$ on the interval $[0,T_{0})$ for some $T_{0}\leq T^{*}$ and we choose $\theta_{M(t)}$ and $\theta_{m(t)}$ such that $d_{\eta}(t)=\eta(\theta_{M(t)},t)-\eta(\theta_{m(t)},t)$ for $t\in[0,T_{0})$. Then, by Lemma 4.2, we get $M(t)\equiv M$ and $m(t)\equiv m$ for some $M,m\in\overline{\mathcal{S}}_{t}$ and thus $d_{\eta}(t)\in\mathcal{C}^{2}((0,T_{0}))$. Thus we obtain from (4.1) that

[TABLE]

for $0<t<T_{0}$, which implies

[TABLE]

due to $|d_{\eta}^{\prime}(0)|\leq d_{v}(0)$. We then define a set

[TABLE]

Due to the assumption on the initial condition, $\mathcal{T}_{1}$ is nonempty. We now claim that $T_{1}^{*}:=\sup\mathcal{T}_{1}=T^{*}$. Suppose, contrary to our claim, that $T_{1}^{*}<T^{*}$. Then, the definition of $T_{1}^{*}$ gives

[TABLE]

The relation (4.6), however, yields

[TABLE]

which is contradictory to (4.7). Thus, we have $T_{1}^{*}=T^{*}$ and the conclusion readily follows. ∎

Using the above uniform boundedenss of the phase diameter, we show the exponential decay of the phase diameter function on the time interval $[0,T^{*})$.

Proposition 4.2**.**

(Exponential decay of phase diameter) Let $T>0$ and $(\eta,h,v)$ be a solution to the system (4.1) on the interval $[0,T]$. Suppose that the initial data satisfy $0<C_{0}<\pi$.

• (i)

If $0<4mkD_{0}<1$, then we have

[TABLE]

where

[TABLE]

• (ii)

If $1\leq 4mkD_{0}$, then we have

[TABLE]

Here $T^{*}>0$ appeared in Lemma 4.2.

Proof.

Similarly as before, we choose $\theta_{M}\in\overline{\mathcal{S}}_{t}$ and $\theta_{m}\in\overline{\mathcal{S}}_{t}$ such that $d_{\eta}(t)=\eta(\theta_{M},t)-\eta(\theta_{m},t)$ on $[0,T^{*})$. Then, by definition, we have

[TABLE]

Then we find from (4.1) that $d_{\eta}$ satisfies

[TABLE]

for $t\in(0,T^{*})$, where we used Lemma 4.3 and

[TABLE]

We now use the relation

[TABLE]

to find

[TABLE]

Using Lemma 2.2 (i), we obtain

[TABLE]

due to $\nu_{1}>\nu_{2}$ and $|d_{\eta}^{\prime}(0)|\leq d_{v}(0)$. This proves $(i)$. The inequality (ii) can also be obtained by applying Lemma 2.2 (ii) to (4.8). This concludes the desired results. ∎

Proposition 4.3**.**

(Exponential decay of velocity diameter) Suppose that $0<C_{0}<\pi$. Then, the following assertions hold.

• (i)

If $0<4mkD_{0}<1$, then we have

[TABLE]

for $t\in[0,T^{*})$, where $C_{1}>0$ is given by

[TABLE]

• (ii)

If $1\leq 4mkD_{0}$, then we have

[TABLE]

for $t\in[0,T^{*})$, where $C_{2}>0$ is given by

[TABLE]

Proof.

In a similar fashion as before, we choose $\theta_{\widetilde{M}}\in\overline{\mathcal{S}}_{t}$ and $\theta_{\widetilde{m}}\in\overline{\mathcal{S}}_{t}$ such that $d_{v}(t)=v(\theta_{\widetilde{M}},t)-v(\theta_{\widetilde{m}},t)$ on the time interval $[0,T^{*})$. Then it follows from (4.1) that $d_{v}$ satisfies

[TABLE]

for $t\in(0,T^{*})$. Applying Grönwall’s lemma to the inequality above gives

[TABLE]

We then now use the upper bounds for $d_{\eta}$ obtained in Proposition 4.2 to have

[TABLE]

for $0<4mkD_{0}<1$, where

[TABLE]

For $1\leq 4mkD_{0}$, we find

[TABLE]

where

[TABLE]

∎

We are now in a position to state the exponential synchronization estimates for the system (4.1) under an appropriate regularity assumptions on the solutions and smallness assumptions on the initial phase and velocity diameters.

Theorem 4.2**.**

Let $T>0$ and $(\eta,h,v)$ be a solution to the system (4.1) on the interval $[0,T]$ with initial data $(h_{0},v_{0})\in H^{2}(\mathcal{S}_{0})\times H^{3}(\mathcal{S}_{0})$. Suppose that there exists $M>0$ independent of $t$ such that

[TABLE]

If the initial phase and velocity diameters $d_{\eta}(0)$, $d_{v}(0)$ are small enough, then there exist positive constants $c_{1},c_{2},\Lambda_{1},\Lambda_{2}$, which are independent of $t$, such that

[TABLE]

for $t\in[0,T)$.

Remark 4.3**.**

Theorem 4.2 implies that as long as there exists a solution satisfying a certain regularity, the exponential decay of phase and velocity diameters can be obtained under some smallness assumptions on the initial phase and velocity diameters. It is worth noticing that we only require the smallness assumptions on the initial phase and velocity diameters, not initial data $(h_{0},v_{0})=(\rho_{0},u_{0})$.

Proof of Theorem 4.2.

It suffices to prove that $T^{*}=T$ in view of Lemma 4.2, Propositions 4.2 and 4.3. Assume to the contrary that $T^{*}<T$. Then, we have

[TABLE]

Recall the notation $v_{c}(t)=\int_{\mathbb{T}}v(\theta,t)\rho_{0}(\theta)\,d\theta$, and note that

[TABLE]

for $t\in[0,T^{*})$. Using the Sobolev embedding, we get

[TABLE]

and furthermore, by Gagliardo-Nirenberg interpolation inequality, we estimate

[TABLE]

for $t\in[0,T^{*})$. This together with (4.9) and (4.11) asserts

[TABLE]

Note that it is clear from the smallness assumption on $d_{\eta}(0)$ and $d_{v}(0)$, and the estimates for $d_{v}(t)$ in Proposition 4.3 that $d_{v}^{2/3}$ and $d_{v}^{2/5}$ are integrable in $(0,\infty)$. More precisely, we obtain

[TABLE]

where $C(d_{\eta}(0),d_{v}(0))>0$ is independent of $t$, and it satisfies

[TABLE]

Thus for sufficiently small $d_{\eta}(0)$ and $d_{v}(0)$, we have

[TABLE]

This is a contradiction to (4.10) and completes the proof. ∎

Remark 4.4**.**

Theorem 4.2 implies

[TABLE]

for some positive constants $c_{3}$ and $\Lambda_{3}$, where $d_{BL}$ denotes the bounded Lipschitz distance. Indeed, if we set

[TABLE]

and $\eta_{\infty}:=\eta_{c}(0)+mv_{c}(0)$, then we have

[TABLE]

for $\phi\in\mathbb{D}$, where $\eta_{c}$ is the mean phase given by

[TABLE]

The last integral can be estimated as follows. First, we get

[TABLE]

For the second term, we use Lemma 4.1 (i) and the relation $\eta_{c}^{\prime}(t)=v_{c}(t)$ to obtain

[TABLE]

Thus we have

[TABLE]

Finally, we use Theorem 4.2 to conclude the desired result.

5. Critical thresholds phenomena

In this section, we study critical thresholds phenomena in the system (1.1). We first provide the local-in-time existence and uniqueness of solutions to the system (1.1).

Theorem 5.1**.**

For any $0<N<M$, there is a positive $T_{0}>0$ such that if

[TABLE]

then the Cauchy problem (1.1)-(1.2) has a unique strong solution

[TABLE]

satisfying

[TABLE]

Proof.

We notice that the local-in-time existence theory is well developed by now, however our solution space is weighted by the distribution function $g$ for natural frequencies. Compared to Theorem 3.1, we need to be more careful because of the convection term in (1.1) which is nonlinear, and it does not appear in the Lagrangian system (1.7). For these reasons, we provide some details of the proof in Appendix A. ∎

Differentiating the momentum equation in (1.1) with respect to $\theta$ and letting $d:=\partial_{\theta}u$, we rewrite the equation (1.1) as follows:

[TABLE]

where $D_{t}$ denotes the time derivative along the characteristic flow $\eta(\theta,\Omega,t)$, i.e., $D_{t}=\partial_{t}+u\partial_{\theta}$.

Proposition 5.1**.**

Consider the system (5.2). Then the following assertions hold.

• (i)

(Subcritical region)* If $1\geq 4Km$ and*

[TABLE]

then $d(\eta(\theta,\Omega,t),\Omega,t)$ remains bounded from below for $(\theta,\Omega)\in\mathbb{T}\times\mathbb{R}$ and $t\geq 0$.

• (ii)

(Supercritical region)* If*

[TABLE]

then $d(\eta(\theta,\Omega,t),\Omega,t)\to-\infty$ in finite time.

Proof.

(i)  Note that the interaction term in (5.2) is easily bounded by

[TABLE]

due to Lemma 2.1 (i), (1.13), and (1.14). This yields

[TABLE]

Let us consider the first inequality in (5.3). It can be rewritten as

[TABLE]

where $d_{\pm}$ is given by

[TABLE]

For $d_{-}\leq d_{0}<d_{+}$, $d(t)\geq d_{0}$ by continuity argument. We now let $q$ solve the following Riccati’s equation:

[TABLE]

The solution of this equation is explicitly given as follows.

[TABLE]

and it is easy to see that $q(t)\geq d_{+}$ if $d_{0}\geq d_{+}.$ The comparison principle yields $d(t)\geq q(t)\geq d_{+}$ for $d_{0}\geq d_{+}$. Thus, we have

[TABLE]

for $t\geq 0$.

(ii) In a similar fashion as above, for the second inequality in (5.3), we get

[TABLE]

This together with the continuity argument implies that if $d_{0}<d_{-}^{*}$, then $d(t)\leq d_{0}<d_{-}^{*}$. This readily gives

[TABLE]

and, subsequently, solving the above differential inequality yields

[TABLE]

Therefore, $d(t)$ will diverge to $-\infty$ until the time $T<(d_{-}^{*}-d_{0})^{-1}$. This completes the proof. ∎

Remark 5.1**.**

In the subcritical case, it follows from the continuity equation in $\eqref{cri_Ku}$ that

[TABLE]

Thus, $\rho$ cannot attain $+\infty$ in a finite time. On the other hand, for the supercritical case, we see

[TABLE]

and thus $\rho$ can be estimated as

[TABLE]

This implies $\rho$ diverges to $+\infty$ until the time $T<(d_{-}^{*}-d_{0})^{-1}$.

Remark 5.2**.**

In the case of no interactions between oscillators, i.e., $K=0$, the momentum equation in (5.2) reduces to the damped pressureless Euler system:

[TABLE]

Thus we obtain a sharp critical thresholds:

• •

If $d_{0}(\theta,\Omega)<-1/m$, then $d(\eta(\theta,\Omega,t),\Omega,t)\to-\infty$ in finite time.

• •

If $d_{0}(\theta,\Omega)\geq-1/m$, then $d(\eta(\theta,\Omega,t),\Omega,t)$ remains bounded for $(\theta,\Omega)\in\mathbb{T}\times\mathbb{R}$ and $t\geq 0$.

We next provide a priori estimates of solutions $(\rho,u)$ to the system (1.1). In the proposition below, we show that $\|\rho\|_{H^{2}_{g}}$ and $\|u\|_{H^{3}_{g}}$ can be controlled by $\|\rho\|_{L^{\infty}}$ and $\|\partial_{\theta}u\|_{L^{\infty}}$.

Proposition 5.2**.**

Let $s\geq 0$ be an integer and consider the system (1.1). Then, for any $T>0$, we have

[TABLE]

Proof.

Replacing $\bar{u}$ in the proof of Theorem 5.1 by $u$, we obtain from (A.4) and (A.9) that

[TABLE]

and

[TABLE]

On the other hand, we find

[TABLE]

due to (1.13). This together with combining (A.4) and (A.9) gives

[TABLE]

Applying Grönwall’s lemma to the inequality above, we conclude the desired result. ∎

As a direct consequence of Theorem 5.1, Propositions 5.1 and 5.2, we have the following results for the critical thresholds phenomena in (1.1).

Theorem 5.2**.**

Let $T>0$ and consider the system (1.1). Under the assumptions of Theorem 5.1, we have the following assertions.

• (i)

(Subcritical region)* If $1\geq 4Km$ and*

[TABLE]

then the system has a global solution, $(\rho,u)\in\mathcal{C}([0,T);H^{2}_{g}(\mathbb{T}\times\mathbb{R}))\times\mathcal{C}([0,T);H^{3}_{g}(\mathbb{T}\times\mathbb{R}))$.

• (ii)

(Supercritical region)* If there exists $(\theta_{*},\Omega_{*})$ such that*

[TABLE]

then the solution blows up in finite time.

Remark 5.3**.**

For reasons mentioned before, almost the same argument as above can be applied to the case of identical oscillators, i.e., $g(\Omega)=\delta_{\Omega_{0}}(\Omega)$ by replacing the weighted spaces $L^{p}_{g}(\mathbb{T}\times\mathbb{R})$ and $H^{s}_{g}(\mathbb{T}\times\mathbb{R})$ by $L^{p}(\mathbb{T})$ and $H^{s}(\mathbb{T})$, respectively.

Remark 5.4**.**

The results of Theorem 5.2 (ii) and Remark 5.1 (ii) give some possible finite-time synchronization. It is very hard to expect the finite-time synchronization phenomena in the classical Kuramoto models, for instances (1.3) and (1.4), with a smooth distribution function $g(\Omega)$ for natural frequencies. However, as mentioned in Introduction, our hydrodynamic model (1.1) is the pressureless Euler-type system, and thus it may form singularities in finite time. It is unclear though how to rigorously justify that this finite-time blow-up of solutions implies the finite-time synchronization. With regard to this matter, we will numerically examine the time evolution of solutions to the system (1.1) in the next section.

6. Numerical experiments

In this section, we present several numerical experiments validating our theoretical results for the system (1.1). We also numerically examine that our system (1.1) exhibits the phase transitions and hysteresis phenomena like the particle system (1.3). For the numerical integration of the system (1.1), the finite volume method is used, in particular, we employ Kurganov-Tadmor central scheme proposed in [28] for the evaluation of numerical fluxes. A brief description of the scheme is provided below.

6.1. Numerical scheme

Note that the system (1.1) can be written in the following form.

[TABLE]

Here we set

[TABLE]

and the source term is given by

[TABLE]

The cell average $Q_{j}$ over the grid cell $\mathcal{C}_{j}=(\theta_{j-1/2},\theta_{j+1/2})$ at time $t$ is given by

[TABLE]

where $\Delta\theta=\theta_{j+1/2}-\theta_{j-1/2}$. Then, for given $\Omega$, the finite volume method is formulated as follows:

[TABLE]

Here, $F^{*}_{j\pm 1/2}$ denotes the numerical flux through the cell interface at $\theta_{j\pm 1/2}$, which will be given later, and $G_{j}^{*}$ is an associated source term evaluated at $\theta=\theta_{j}$ where the integration is performed using the midpoint rule. The reconstruction first requires a piecewise linear function

[TABLE]

where $\theta_{j}$ is the center of $\mathcal{C}_{j}$, and $\sigma_{j}^{n}$ denotes an approximation to the spatial derivative on $\mathcal{C}_{j}$. In order to prevent nonphysical oscillations, we use the slope limiter method which was introduced by van Leer [39]. In particular, we use minmod slope here:

[TABLE]

where the minmod function is defined by

[TABLE]

The numerical fluxes are now evaluated as

[TABLE]

Here, $a^{+}_{j+1/2}$ and $a^{-}_{j+1/2}$ denote the largest and the smallest speed of characteristic at the cell interfaces. The reconstructed values at the cell interface $\theta_{j+1/2}$ using (6.1) are given by

[TABLE]

Finally, the second-order Runge–Kutta method is employed for time integrations.

6.2. Time evolutions of density and velocity

In this subsection, we present the time evolutions of the density and the velocity profiles at different time $t$’s for both identical and nonidentical oscillators case. For the identical case, the $\Omega$-dependences of $\rho$ and $u$ are simply neglected, see Section 4. We identify $\mathbb{T}$, the $\theta$-domain, as $[-\pi,\pi]$ for the numerical computation domain and set the initial density $\rho_{0}$ and the distribution function $g(\Omega)$ the standard normal distribution. The numbers of $\theta$-grid and $\Omega$-grid are $1000$ and $600$, respectively.

Figure 1 exhibits the time evolutions of $\rho$ and $u$ for the identical cases. The parameters and the initial data for Figures 1 (a) and (b) are selected as

[TABLE]

so that they lie in subcritical region. We can observe the contractivity of support of $\rho$, and it is consistent with our theoretical result, Remark 4.4. Nevertheless, $\rho$ does not blow-up in finite time, see Remark 5.1. We also note that the velocity diameter $d_{v}$ decreases. On the other hand, Figures 1 (c) and (d) show the profiles of $\rho$ and $u$ for supercritical region around $\theta=0$ with the following initial velocity and $m,K$:

[TABLE]

As proved in Proposition 5.1 and Remark 5.1, we can observe the finite-time blow-up of solutions. Furthermore, if we set

[TABLE]

so that we have a supercritical case around $\theta=\pm\pi/2$, the profiles show the finite-time blow-up around these points as can be seen in Figures 1 (e) and (f). It is remarkable that the density forms Dirac measures in finite time in the supercritical case, and this supports that our model exhibits the finite-time synchronization phenomena.

We next consider the nonidentical case, that is, the distribution function $g$ for natural frequencies is not the form of Dirac measure centered on some fixed point $\Omega_{0}\in\mathbb{R}$.

Figures 2 (a), (b) and Figure 3 illustrate the profiles for the subcritical cases in the two-dimensional and three-dimensional plots, respectively. In Figures 2 (a) and (b), we plot modified density $\tilde{\rho}$ and velocity $\tilde{u}$ which are given by

[TABLE]

so that the profiles can be observed on $\mathbb{T}$. The parameters and the initial velocity are set

[TABLE]

which guarantees the subcritical region for all $(\theta,\Omega)$ so that we have the global existence of the solution, see Theorem 5.2 (i). We can observe that the density does not tend to a Dirac measure unlike identical case. Note that for each fixed $\Omega$, the density profile $\rho(\theta,\Omega)$ exhibits the bell-shaped distribution whose center varies according to the natural frequency $\Omega$. The frequency synchronization is also observed in the time evolution of $u$.

The plots of $\tilde{\rho},\rho,\tilde{u},u$ for the first supercritical case are provided in Figures 2 (c), (d) and Figure 4 where the parameters and the initial data are chosen as

[TABLE]

In this case, the initial data lie in the supercritical region around $\theta=0$. The finite-time blow-up of $\rho$ and $\partial_{\theta}u$ in the small time interval is easily observed here, which is consistent with our theoretical results Proposition 5.1 and Remark 5.1.

In Figures 2 (e), (f) and Figure 5, we present the plots of the density and velocity profiles for the second supercritical case. As is the case in the identical oscillators(Figures 1 (e) and (f)), we take the following parameters and the initial velocity:

[TABLE]

to consider the supercritical case around $\theta=\pm\pi/2$, and the figures also exhibit the finite-time blow-up around these points.

6.3. Phase transitions & hysteresis phenomena

In Figure 6, we show the phase transitions of the order parameter $r^{\infty}$ for the hydrodynamic model (1.1) with the coupling strength on the interval $[0,4]$. It is known that unlike the Kuramoto model without inertia, where the phase transition of the order parameter versus the coupling strength is continuous provided the distribution function $g$ is Gaussian, even the small inertia can lead to a discontinuous and hysteretic phase transition for the system (1.3). As noted in Introduction, the hydrodynamic model (1.1) also carries the hysteresis phenomena as (1.3) does. Figure 6 exhibits the discontinuous phase transition of (1.1) with $m=0.1,0.5,1$. The numbers of $\theta$-grid and $\Omega$-grid are set to $100$ and $600$, respectively. We increase $K$ from [math] to $4$ with the mesh spacing of $0.1$ for $K$. When $K$ reaches $4$, the same procedure is iterated by decreasing $K$ back to [math]. In order to gain the clear observation on the thresholds, the finer mesh spacing in $K$, which is $0.05$, is used around them. The direction of jump is indicated with arrows. The initial density $\rho_{0}$ and the distribution $g$ are set to the Gaussian distribution and the initial velocity is $u_{0}(\theta,\Omega)=-(0.5)\sin\theta$. Note that $u_{0}$ is chosen such that it does not lie in the supercritical region for any of $m\in\{0.1,0.5,1\}$ and $K\in[0,4]$, see Section 5. We can observe that as $m$ increases, the hysteresis becomes more noticeable.

7. Conclusion

In this manuscript, we presented a new hydrodynamic model for the synchronization phenomena and discussed the local-in-time existence theory. For the identical natural frequencies, we provided two different approaches for the synchronization estimates; kinetic energy combined with the order parameter estimates and the second-order Grönwall-type inequality estimates on the phase and velocity diameters. In particular, by the latter strategy, we showed that the limiting density is the form of the Dirac measure. We also analyzed the critical threshold phenomena in our main system. By this analysis, we found that classical solutions can be blow-up in finite time, which is not observed in the classical Kuramoto models. We were not able to prove this finite-time blow-up of solutions implies the finite-time synchronization, however, numerical simulations illustrated that the density with initial data in the supercritical region becomes Dirac measures in finite time. We also presented several numerical simulations validating our analytical results. The numerical results showed that our main system also has similar features, such as phase transitions and hysteresis phenomena, compared to the Kuramoto model with inertia. As briefly mentioned in Introduction, the pressureless Euler-type equations may develop the formation of singularities. For this reason, it is natural to take into account the notion of measure-valued solutions. Thus it would be interesting to study the existence of measure-valued solutions to our main system. This may enable us to have the global-in-time regularity of solutions. As the first step in this hydrodynamic modeling of synchronization phenomena, we only deal with the case of identical oscillators for the synchronization estimates. Hence, our next step is to generalize our analysis for the case of nonidentical natural frequencies. We will investigate these interesting issues in future.

Appendix A Proof of Theorem 5.1

For computational simplicity, we set $m=K=1$. Let $T>0$ be given, and we consider the system:

[TABLE]

with the initial data $(\rho_{0},u_{0})$ satisfying the assumptions (5.1). Here $\bar{u}$ satisfies

[TABLE]

Note that we can use a standard linear theory to show the existence of solutions to the system (A.1). We begin by estimating $\|\rho\|_{H^{2}_{g}}$. A direct calculation gives

[TABLE]

Similarly, we can easily obtain

[TABLE]

We next estimate

[TABLE]

Note that

[TABLE]

then this together with applying Hölder’s inequality yields

[TABLE]

Thus we get

[TABLE]

Using the above inequality, we estimate the last term on the right hand side of (A.3) as

[TABLE]

We now combine all of the above observations to have

[TABLE]

that is,

[TABLE]

On the other hand, by taking into account the characteristic flow defined by

[TABLE]

we can easily estimate

[TABLE]

for $t\in[0,T]$. This together with (A.4) asserts

[TABLE]

Applying Grönwall’s lemma to the inequality above, we have

[TABLE]

for $t\in[0,T]$.

We next estimate $\|\partial_{\theta}u\|_{L^{\infty}}$ and $\|u\|_{H^{3}_{g}}$. For this, similarly as before, we can use the characteristic flow to find

[TABLE]

for all $(\theta,\Omega)\in\mathbb{T}\times\mathbb{R}$. This enables us to divide the momentum equation in (A.1) by $\rho$, and this gives

[TABLE]

Along the characteristic flow (A.5), we estimate

[TABLE]

Thus we obtain

[TABLE]

for $t\in[0,T]$, where $C>0$ is independent of $t$. For the estimate of $\|u\|_{H^{3}_{g}}$, we first notice that

[TABLE]

for $k\in\mathbb{N}$, due to Remark 2.1. This together with similar estimates as above yields

[TABLE]

where we used the assumption $\int_{\mathbb{R}}\Omega^{2}g(\Omega)\,d\Omega<\infty$ in (1.13) and the following inequality:

[TABLE]

This asserts

[TABLE]

Here we used the estimate (A.8) and the assumption (A.2). Applying Grönwall’s lemma to (A.9) gives

[TABLE]

By combining this with (A.6), (A.7), and (A.8), we have

[TABLE]

We finally choose $T>0$ small enough such that the right hand side of the above inequality is less than $M$. We then deal with the approximations for the system (1.1):

[TABLE]

with the initial data and the first iteration step given by

[TABLE]

and

[TABLE]

For the system (A.10), we use the previous argument to get

[TABLE]

for some $T_{0}>0$. We next show that $(\rho^{n},u^{n})$ is a Cauchy sequence in $\mathcal{C}([0,T_{0}];L^{2}_{g}(\mathbb{T}\times\mathbb{R}))\times\mathcal{C}([0,T_{0}];H^{1}_{g}(\mathbb{T}\times\mathbb{R}))$. For this, we introduce the following simplified notations:

[TABLE]

Then straightforward computations yield

[TABLE]

where we used

[TABLE]

due to (A.11). This asserts

[TABLE]

We next estimate $\|u^{n+1,n}\|_{L^{4}_{g}}$ and $\|\partial_{\theta}u^{n+1,n}\|_{L^{2}_{g}}$. It follows from the momentum equation in (A.10) that

[TABLE]

Thus we find

[TABLE]

where we used (1.13). This gives

[TABLE]

We also use the similar argument as above to estimate

[TABLE]

Here we used

[TABLE]

Thus we obtain

[TABLE]

and this together with (A.12) and (A.13) yields

[TABLE]

for $t\in[0,T_{0}]$, where $C>0$ is independent of $n$. Since

[TABLE]

this conclude that $(\rho^{n},u^{n})$ is a Cauchy sequence in $\mathcal{C}([0,T_{0}];L^{2}_{g}(\mathbb{T}\times\mathbb{R}))\times\mathcal{C}([0,T_{0}];H^{1}_{g}(\mathbb{T}\times\mathbb{R}))$. The rest part of the proof is almost the same with Steps 3 - 5 in the proof of Theorem 3.1. This completes the proof.

Acknowledgments

YPC was supported by National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIP) (No. 2017R1C1B2012918) and POSCO Science Fellowship of POSCO TJ Park Foundation.