# Mean field Kuramoto models on graphs

Wuchen Li
Department of Mathematics
University of South Carolina, Columbia, SC 29208, USA<EMAIL_ADDRESS>and  Hansol Park
Department of Mathematics
Simon Fraser University, 8888 University Dr, Burnaby, BC V5A 1S6, Canada
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

One of a classical synchronization model is the Kuramoto model. We propose
both first and second order Kuramoto dynamical models on graphs using discrete
optimal transport dynamics. We analyze the synchronization behaviors for some
examples of Kuramoto models on graphs. We also provide a generalized Hopf-Cole
transformation for discrete optimal transport systems. Focus on the two points
graph, we derive analytical formulas of the Kuramoto dynamics with various
potential induced from entropy functionals. Several numerical examples for the
Kuramoto model on general graphs are presented.

###### Key words and phrases:

Kuramoto model; Emergent; Optimal transport; Schrödinger bridge problem;
Graph; Hopf-Cole transformation on graphs.

###### 2020 Mathematics Subject Classification:

34D06, 49Q22

Acknowledgment. W. Li is supported by a start up funding in University of
South Carolina and NSF RTG: 2038080. And H. Park is supported by Pacific
Institute for the Mathematical Science(PIMS), Canada postdoctoral fellowship.

##  1\. Introduction

Synchronization phenomena of interacting particles plays essential roles in
physics, biology, social systems, and artificial intelligence (AI). This field
of research have been studied intensively by two pioneers: A. Winfree [21] and
Y. Kuramoto [12, 13]. One of the famous synchronization model is the Kuramoto
model. This model is a synchronization model on the unit circle given as
follows:

(1.1)
$\displaystyle\begin{cases}\displaystyle\dot{\varphi}_{i}=\nu_{i}+\frac{\kappa}{N}\sum_{k=1}^{N}\sin(\varphi_{k}-\varphi_{i})\quad\forall~{}t>0,\\\
\varphi_{i}(0)=\varphi_{i}^{0}\in\mathbb{R}\quad\forall~{}i=1,2,\dots,N.\end{cases}$

The emergent behavior of system (1.1) has been studied in [9, 11]. For the
identical case (i.e. $\nu_{i}=0$ for all $i$), we can prove that the
difference between phases $\varphi_{i}-\varphi_{j}$ converges to $2\pi m$
exponentially for some $m\in\mathbb{Z}$ under suitable initial conditions. If
the number of particles $N$ tends to infinity, we have the following mean
field dynamical model of the identical case of (1.1):

(1.2)
$\displaystyle\begin{cases}\partial_{t}\rho_{t}(\varphi)+\partial_{\varphi}(\rho_{t}(\varphi)L[\rho_{t}](\varphi))=0,\\\
\displaystyle
L[\rho_{t}](\varphi)=-\kappa\int_{\mathbb{S}^{1}}\sin(\varphi-\varphi_{*})\rho_{t}(\varphi_{*})d\varphi_{*},\quad\varphi\in\mathbb{R},\end{cases}$

where $\rho_{0}$ is the initial distribution. Similarly, the emergent behavior
of equation (1.2) is given as

(1.3) $\displaystyle\rho_{t}\to\delta_{\varphi^{\infty}}$

in [10], where $\delta_{\varphi^{\infty}}$ is a Dirac delta function
concentrated at a constant point $\varphi^{\infty}$. We can interpret it as
the probability density function defined on the unit circle is concentrating
into one point $\varphi^{\infty}$ as time goes to infinity. In practice, a
large population of particles are often settled on a discrete spatial domain;
e.g. a simple finite graph. A natural question arise:

(Q) What is the synchronization phenomena of a large population of interacting
particles on a discrete domain, such as finite graphs?

This paper provides a proper form of a discretized system (1.2) on finite
graphs using optimal transport dynamics [20, 8, 18]. We show that if we
discretize system (1.2) onto a discrete graph with $n$ vertices, then the
corresponding emergent behavior to (1.3) can be expressed as follows (up to
permutation):

$\displaystyle(\rho_{1}(t),\rho_{2}(t)\cdots,\rho_{n}(t))\to(1,0,\cdots,0)\quad\text{as}\quad
t\to\infty.$

To study discrete dynamical system (1.2), we use the tools of optimal
transport dynamics on graphs; see [14, 15, 17, 7]. In Section 2.2, we express
system (1.2) as both a gradient flow system and a Hamiltonian system. From
this fact, we have the following first order model on the discrete model as
follows:

$\begin{cases}\displaystyle\frac{d\rho_{j}}{dt}=\kappa\sum_{k=1}^{n}\theta_{jk}(\rho_{j}-\rho_{k}),\quad
t>0,\\\ \rho_{j}(0)=\rho_{j}^{0},\quad j=1,2,\cdots,N.\end{cases}$

Also, we can follow the discretization introduced in [5, 6] to obtain the
discretized system of (1.2) below:

$\displaystyle\begin{cases}\displaystyle\frac{d\rho_{j}}{dt}-\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-S_{l})\theta_{jl}=0,\vspace{0.2cm}\\\
\displaystyle\frac{dS_{j}}{dt}+\frac{1}{2}\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-S_{l})^{2}\frac{\partial\theta_{jl}}{\partial\rho_{j}}=\frac{\partial}{\partial\rho_{j}}\left(\frac{1}{2}\sum_{(k,l)\in
E}\frac{\omega_{kl}}{2}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})^{2}\theta_{kl}\right),\vspace{0.2cm}\\\
\displaystyle\rho_{j}(0)=\rho_{j}^{0},\quad S_{j}(0)=S_{j}^{0},\quad
S_{j}(0)=-\frac{\partial\mathcal{F}(\rho)}{\partial\rho_{j}}\Bigg{|}_{\rho=\rho(0)},\end{cases}$

where $\mathcal{F}(\rho)$ is a potential functional defined on a graph,
$\rho_{j}$ is the probability density on node $j$, and $S_{j}$ is the
potential function on node $j$. Also, $\theta_{ij}$ represents the probability
weight of the edge connecting $i^{th}$ and $j^{th}$ vertices. For example,
$\theta_{ij}=\frac{1}{2}(\rho_{i}+\rho_{j})$ can be picked, however, in this
paper we consider the general probability weight given as
$\theta_{ij}=\theta(\rho_{i},\rho_{j})$ for some function $\theta$. We study
this discretized model in this paper intensively.

The main results of this paper are three-fold. First, we develop first and
second order Kuramoto model on graphs following discrete optimal transport
dynamics. We also provide emergent behaviors (synchronization behaviors) of
the Kuramoto model on graphs. Second, we provide the Hopf-Cole transformation
for discrete optimal transport system. Third, we derive analytical formulas of
the Kuramoto dynamics on the two points graph with various entropy
functionals.

This paper is organized as follows. In Section 2, we review some aggregation
models, optimal transport dynamics on discrete graph, and the generalized
Hopf-Cole transformation. We construct first and second order Kuramoto model
on finite graphs and study their emergent behaviors in Section 3. We also
apply the generalized Hopf-Cole transformation on discrete graph in Subsection
3.4. In Section 4, we provide analytical formulas of the Kuramoto dynamics on
the two points graph. Several numerical examples are presented in Section 5.

##  2\. Review

In this section, we provide some reviews on concepts and results on
aggregation models and optimal transport equations.

For the intuition, consider a case of classical mechanics on $\mathbb{R}^{d}$.
We first consider a gradient flow system

$\begin{cases}\dot{x}=-\nabla U(x),\\\
x(0)=x^{0}\in\mathbb{R}^{d}.\end{cases}$

We second study a Hamiltonian system. Suppose the Hamiltonian is

$H(x,p)=\frac{1}{2}p^{2}-\frac{1}{2}|\nabla U(x)|^{2}.$

The Hamiltonian system satisfies

(2.1)
$\displaystyle\begin{cases}\displaystyle\dot{x}=p,\quad\dot{p}=-\frac{1}{2}\nabla|\nabla
U(x)|^{2},\\\ x(0)=x^{0},\quad p(0)=p^{0}.\end{cases}$

If we impose the following assumption:

$p^{0}=-\nabla U(x^{0}),$

then this Hamiltonian system can be reduced to $\dot{x}=-\nabla U(x)$. See
Lemma A.1 for detail proof. In this sense, any gradient flow system can be
expressed as a particular case of Hamiltonian system.

###  2.1. Review on aggregation models

In this subsection, we review some gradient flow type aggregation models. The
Kuramoto model is the simplest synchronization model defined on
$\mathbb{S}^{1}$. Consider

(2.2)
$\displaystyle\begin{cases}\displaystyle\dot{\varphi}_{i}=\nu_{i}+\frac{\kappa}{N}\sum_{k=1}^{N}\sin(\varphi_{k}-\varphi_{i}),\\\
\varphi_{i}(0)=\varphi_{i}^{0},\quad\forall~{}i\in\mathcal{N}:=\\{1,2,\cdots,N\\},\end{cases}$

where $N$ is the number of particles, $\nu_{i}$ is the natural frequency of
the $i^{th}$ particle, $\kappa$ is the coupling strength, and $\varphi_{i}$ is
the phase of the $i^{th}$ particle. If the natural frequency is identical,
i.e., $\nu_{i}\equiv\nu$ for all $i\in\mathcal{N}$, then solution of (2.2)
$\Phi:=\\{\varphi_{i}\\}_{i=1}^{N}$ exhibits the complete synchronization:

$\lim_{t\to\infty}(e^{\mathrm{i}\varphi_{i}}-e^{\mathrm{i}\varphi_{j}})=0.$

Also, if the coupling strength $\kappa$ is large enough, then

$\lim_{t\to\infty}(\dot{\varphi}_{i}-\dot{\varphi}_{j})=0.$

The Kuramoto model can be expressed as a gradient flow system:

$\displaystyle\begin{cases}\dot{\Phi}=-\nabla_{\Phi}\mathcal{V},\quad t>0,\\\
\displaystyle\mathcal{V}(\Phi)=-\sum_{i=1}^{N}\nu_{i}\varphi_{i}-\frac{\kappa}{2N}\sum_{i,j=1}^{N}\cos(\varphi_{i}-\varphi_{j}).\end{cases}$

When $N$ goes to infinity, we can obtain the mean field version of system
(2.2). The mean field version of the identical Kuramoto model($\nu_{i}\equiv
0$) is the Kuramoto-Sakaguchi model [19] given as follows:

(2.3)
$\displaystyle\begin{cases}\partial_{t}\rho_{t}(\varphi)+\partial_{\theta}(\rho_{t}(\varphi)L[\rho_{t}](\varphi))=0,\\\
\displaystyle
L[\rho_{t}](\varphi)=-\kappa\int_{\mathbb{S}^{1}}\sin(\varphi-\varphi_{*})\rho_{t}(\varphi_{*})d\varphi_{*}\quad\forall~{}\varphi\in\mathbb{R},\end{cases}$

where $\rho_{0}$ is an initial probability distribution function defined on
$\mathbb{S}^{1}$.

A generalization of system (2.2) to general sphere $\mathbb{S}^{d-1}$ with
$d\geq 2$ is introduced in [16]:

$\begin{cases}\displaystyle\dot{x}_{i}=\Omega_{i}x_{i}+\frac{\kappa}{N}\sum_{k=1}^{N}\left(x_{k}-\langle
x_{i},x_{k}\rangle x_{i}\right),\\\
x_{i}(0)=x_{i}(0)\in\mathbb{S}^{d-1}\subset\mathbb{R}^{d},\quad\forall~{}i\in\mathcal{N},\end{cases}$

where $\Omega_{i}$ is a skew-symmetric matrix of size $d\times d$. We call
this system the swarm sphere model. The mean field version of the identical
swarm sphere model ($\Omega_{i}\equiv 0$) is given below:

(2.4)
$\displaystyle\begin{cases}\partial_{t}\rho_{t}(x)+\nabla_{\mathbb{S}^{d}}\cdot(\rho_{t}(x)L[\rho_{t}](x))=0,\\\
\displaystyle L[\rho](x)=\kappa\int_{\mathbb{S}^{d}}(y-\langle x,y\rangle
x)\rho(y)dy\quad\forall~{}x\in\mathbb{S}^{d},\end{cases}$

where $\rho_{0}$ is the initial distribution on $\mathbb{S}^{d}$. We will use
this system (2.4) to construct the consensus model on a regular $n$-simplex.

###  2.2. Optimal transport dynamics.

In this subsection, we present some previous results on optimal transport
dynamics. We first consider a gradient flow system on Wasserstein space. We
then demonstrate that any Wasserstein gradient flows can be expressed as
Hamiltonian flows.

Consider the following Wasserstein gradient flow system:

(2.5)
$\displaystyle\partial_{t}\rho_{t}=-\mathrm{grad}_{W}\mathcal{F}(\rho_{t})=\mathrm{div}\left(\rho_{t}\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)\right),$

where the initial distribution is $\rho_{0}$, $\mathcal{F}$ is a given
functional, $\frac{\delta}{\delta\rho}$ is the $L^{2}$ first variation
operator w.r.t. $\rho$, and $\mathrm{grad}_{W}$ is a gradient operator in
Wasserstein space:

$\mathrm{grad}_{W}\mathcal{F}[\rho]=-\mathrm{div}\left(\rho\nabla\frac{\delta\mathcal{F}[\rho]}{\delta\rho}\right).$

Consider the Wasserstein Hamiltonian system

(2.6)
$\displaystyle\begin{cases}\partial_{t}\rho_{t}=\partial_{S}\mathcal{H}(\rho_{t},S_{t}),\\\
\partial_{t}S_{t}=-\partial_{\rho}\mathcal{H}(\rho_{t},S_{t}),\end{cases}$

where the Hamiltonian is defined as

$\mathcal{H}(\rho,S)=\frac{1}{2}\int\Big{(}|\nabla
S(x)|^{2}-|\nabla\delta\mathcal{F}(\rho)|^{2}\Big{)}\rho(x)dx.$

In other words,

(2.7)
$\displaystyle\begin{cases}\partial_{t}\rho_{t}+\mathrm{div}(\rho_{t}\nabla
S_{t})=0,\\\ \displaystyle\partial_{t}S_{t}+\frac{1}{2}|\nabla
S_{t}|^{2}=\frac{\delta}{\delta\rho_{t}}\left(\frac{1}{2}\int\left|\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)\right|^{2}\rho_{t}(x)dx\right).\end{cases}$

If the initial data $(\rho_{0},S_{0})$ satisfies

(2.8) $\displaystyle
S_{0}=-\frac{\delta\mathcal{F}(\rho_{0})}{\delta\rho_{0}},$

then system (2.7) can be reduced to (2.5); see Lemma A.2 for the detail proof.
So, we conclude that any Wasserstein gradient flows can also be expressed as
particular Hamiltonian systems.

###  2.3. Generalized Hopf-Cole transformation

In the remained part of this section, we review generalized Hopf-Cole
transformations [2, 14]. It is to rewrite the Hamiltonian system in term of
entropy functions. Define a pair of functions $(\eta,\eta^{*})$, such that

$\displaystyle\begin{cases}\delta\mathcal{F}(\rho)=\delta\mathcal{F}(\eta)+\delta\mathcal{F}(\eta^{*}),\\\
S=\delta\mathcal{F}(\eta)-\delta\mathcal{F}(\eta^{*}).\end{cases}$

The Hamiltonian system (2.6) in term of $(\eta,\hat{\eta})$ satisfies

$\begin{cases}\partial_{t}\eta_{t}=\sigma(\eta_{t},\eta_{t}^{*})\partial_{\eta^{*}}\mathcal{K}(\eta_{t},\eta_{t}^{*}),\\\
\partial_{t}\eta_{t}^{*}=-\sigma(\eta_{t}^{*},\eta)\partial_{\eta}\mathcal{K}(\eta_{t},\eta_{t}^{*}),\end{cases}$

where $\mathcal{K}(\eta,\eta^{*})=\mathcal{H}(\rho,S)$ and

$\sigma(\eta,\eta^{*})(x,w)=-\frac{1}{2}\iint[\delta^{2}\mathcal{F}(\eta)]^{-1}(x,y)\delta^{2}\mathcal{F}(\rho)(y,z)[\delta^{2}\mathcal{F}(\eta^{*})]^{-1}(z,w)dydz.$

This transformation can be not well-defined when the inverse of
$\delta\mathcal{F}$ is not well-defined. One can also set new variables $\xi$
and $\xi^{*}$, such that

(2.9)
$\displaystyle\xi=\delta\mathcal{F}(\eta),\quad\xi^{*}=\delta\mathcal{F}(\eta^{*}).$

Then, the Hamiltonian system (2.8) in term of $(\xi,\xi^{*})$ satisfies

(2.10)
$\displaystyle\begin{cases}\partial_{t}\xi(x)\displaystyle=\nabla\xi(x)\cdot\nabla\xi^{*}(x)-\int[\delta^{2}\mathcal{F}](x,u)\nabla\cdot(\rho(u)\nabla\xi(u))du,\\\
\partial_{t}\xi^{*}(x)\displaystyle=-\nabla\xi(x)\cdot\nabla\xi^{*}(x)+\int[\delta^{2}\mathcal{F}](x,u)\nabla\cdot(\rho(u)\nabla\xi^{*}(u))du.\end{cases}$

For the detail proof, we refer [2]. In this paper, we perform the discretized
analog of the change of variable (2.10) in Section 3.4.

###  2.4. Aggregation models via optimal transport dynamics

Actually, many kinetic aggregation models can be expressed as Wasserstein
gradient flow system (2.5). In detail, the kinetic Kuramoto model (2.3)
satisfies

$\displaystyle\begin{cases}\displaystyle\partial_{t}\rho_{t}(\theta)=\partial_{\theta}\left(\rho_{t}(\theta)\partial_{\theta}\left(\frac{\delta}{\delta\rho_{t}(\theta)}\mathcal{V}(\rho_{t})\right)\right)=-\mathrm{grad}_{W}(\mathcal{V}(\rho_{t})),\\\
\displaystyle\mathcal{V}(\rho)=-\frac{\kappa}{2}\iint\cos(\theta_{*}-\theta_{**})\rho(\theta_{*})\rho(\theta_{**})d\theta_{*}d\theta_{**}.\end{cases}$

Similarly, the kinetic swarm sphere model (2.4) forms

(2.11)
$\displaystyle\begin{cases}\partial_{t}\rho_{t}(x)=\nabla_{\mathbb{S}^{d}}\cdot\Big{(}\rho_{t}(x)\nabla_{\mathbb{S}^{d}}\frac{\delta}{\delta\rho_{t}}\mathcal{V}[\rho_{t}]\Big{)}=-\mathrm{grad}_{W}\mathcal{V}[\rho_{t}](x),\\\
\displaystyle\mathcal{V}[\rho]=\frac{\kappa}{2}\iint_{(\mathbb{S}^{d})^{2}}\|x-y\|^{2}\rho(x)\rho(y)dxdy.\end{cases}$

Recall the kinetic swarm sphere model is given as (2.4) and its gradient flow
formulation can be expressed as (2.11). Since system (2.11) is a special case
of (2.5), we can express this system as a Hamiltonian flow formulation using
(2.7) as follows:

(2.12)
$\displaystyle\begin{cases}\displaystyle\partial_{t}\rho_{t}+\mathrm{div}(\rho_{t}\nabla
S_{t})=0,\\\ \displaystyle\partial_{t}S_{t}+\frac{1}{2}|\nabla
S_{t}|^{2}=\frac{\delta}{\delta\rho_{t}}\left(\frac{1}{2}\int\left|\nabla\left(\frac{\delta\mathcal{V}(\rho_{t})}{\delta\rho_{t}}\right)\right|^{2}\rho_{t}(x)dx\right),\end{cases}$

where $\mathcal{V}$ is defined in (2.11).

##  3\. Main results: Kuramoto models on graphs

In section 3, we present the main result of this paper.

###  3.1. Modeling of synchronization model on discrete graph

We first connect aggregation models with optimal transport dynamics.

Let the regular $n$-simplex be given on
$\mathbb{S}^{n-2}\subset\mathbb{R}^{n-1}$($n\geq 2$). Let $n$ points
$x_{1},x_{2},\cdots,x_{n}$ satisfy

$\\{x_{1},\cdots,x_{n}\\}\subset\mathbb{S}^{n-2},\quad\|x_{i}-x_{j}\|=L_{n}\quad\forall~{}i\neq
j$

for some positive constant $L_{n}$. Actually, this constant can be calculated
for all $n\geq 2$. e.g., $n=3$, it is an equilateral triangle with side length
$L_{3}=\sqrt{3}$. We consider a graph $G_{n}$, which consists
$\\{x_{1},\cdots,x_{n}\\}$ as the set of vertices and $\\{(x_{i},x_{j}):1\leq
i<j\leq n\\}$ as the set of edges. Since (2.12) is a special case of (2.7), we
have the following discretization on $G_{n}$ of system (2.12):

(3.1)
$\displaystyle\begin{cases}\displaystyle\frac{d\rho_{j}}{dt}-\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{n}\omega_{jl}(S_{j}-S_{l})\theta_{jl}=0,\\\
\displaystyle\frac{dS_{j}}{dt}+\frac{1}{2}\sum_{\begin{subarray}{c}l=1\\\
l\neq
j\end{subarray}}^{n}(S_{j}-S_{l})^{2}\frac{\partial\theta_{jl}}{\partial\rho_{j}}=\frac{\partial}{\partial\rho_{j}}\left(\frac{1}{4}\sum_{\begin{subarray}{c}k,l=1\\\
k\neq
l\end{subarray}}^{n}\omega_{kl}\Big{(}\partial_{k}\tilde{\mathcal{F}}_{n}-\partial_{l}\tilde{\mathcal{F}}_{n}\Big{)}^{2}\theta_{kl}\right),\end{cases}$

where potential $\tilde{\mathcal{F}}$ is given as

$\tilde{\mathcal{F}}_{n}(\rho)=\frac{\kappa}{2}\sum_{i,j=1}^{n}\|x_{i}-x_{j}\|^{2}\rho_{i}\rho_{j}=\frac{\kappa}{2}\sum_{i,j=1}^{n}L_{n}^{2}(1-\delta_{ij})\rho_{i}\rho_{j}.$

Here $\rho_{j}$ and $S_{j}$ are the probability densityh and potential
function on node $j$, respectively. Since the potential function will be
differentiated, we add some constant to define another potential function:

(3.2)
$\displaystyle\mathcal{F}_{n}[\rho]=\tilde{\mathcal{F}}_{n}[\rho]-\frac{\kappa}{2}L_{n}^{2}=-\frac{\kappa}{2}\sum_{i,j=1}^{n}\delta_{ij}\rho_{i}\rho_{j}=-\frac{\kappa
L_{n}^{2}}{2}\sum_{i=1}^{n}\rho_{i}^{2}.$

Also, from $\|x_{i}-x_{j}\|\equiv L_{n}$ for all $i\neq j$, we set
$\omega_{ij}=\frac{1}{L_{n}^{2}}$ for all $i\neq j$. We rescale the system to
set $L_{n}=1$. For the convenience, we set $\theta_{ii}=\rho_{i}$. From the
definition of $\mathcal{F}_{n}$, we get

$\partial_{i}\mathcal{F}_{n}=-\kappa\rho_{i}\quad\forall~{}1\leq i\leq n.$

From these additional assumptions, we can reduce system (3.1) as follows:

(3.3)
$\displaystyle\begin{cases}\displaystyle\frac{d\rho_{j}}{dt}-\sum_{l=1}^{n}(S_{j}-S_{l})\theta_{jl}=0,\\\
\displaystyle\frac{dS_{j}}{dt}+\frac{1}{2}\sum_{l=1}^{n}(S_{j}-S_{l})^{2}\frac{\partial\theta_{jl}}{\partial\rho_{j}}=\frac{\partial}{\partial\rho_{j}}\left(\frac{\kappa^{2}}{4}\sum_{k,l=1}^{n}\Big{(}\rho_{k}-\rho_{l}\Big{)}^{2}\theta_{kl}\right).\end{cases}$

If the initial data of system (3.3) satisfy

(3.4) $\displaystyle
S_{j}(0)=-\partial_{j}\mathcal{F}_{n}\Big{|}_{t=0}=\kappa\rho_{j}(0)\quad\forall~{}1\leq
j\leq n,$

then the system can be reduced to the gradient flow formulation(first order
system).

###  3.2. First order dynamics

In this subsection, we study the emergent behavior of system (3.3) with (3.4).
From Lemma A.2, we have

$S_{j}(t)=\kappa\rho_{j}(t)$

for all $1\leq j\leq n$ and $t\geq 0$. We substitute this relation into (3.3),
then we have

(3.5)
$\displaystyle\frac{d\rho_{j}}{dt}=\kappa\sum_{k=1}^{n}\theta_{jk}(\rho_{j}-\rho_{k}).$

From a simple calculation, we get

$\frac{d}{dt}\sum_{j=1}^{n}\rho_{j}^{2}=\kappa\sum_{j,k=1}^{n}\theta_{jk}\rho_{j}(\rho_{j}-\rho_{k}).$

We assume the following conditions on $\theta_{ij}=\theta(\rho_{i},\rho_{j})$:

$\bullet$ Symmetry: $\theta_{ij}=\theta_{ji}$,

$\bullet$ Lipschitz continuity: $\theta_{ij}$ is Lipschitz continuous
function,

$\bullet$ $\theta_{ij}\geq 0$ and the equality only holds either $\rho_{i}=0$
or $\rho_{j}=0$.

Then, we have

$\frac{d}{dt}\sum_{j=1}^{n}\rho_{j}^{2}=\frac{\kappa}{2}\sum_{j,k=1}^{n}\theta_{jk}(\rho_{j}-\rho_{k})^{2}\geq
0.$

On the other hand, we can easily check that if $\theta_{jk}$ is continuous,
then the derivative of $\frac{d}{dt}\sum_{j=1}^{n}\rho_{j}^{2}$ is uniformly
continuous and we can apply Barbalat’s lemma to get

$\lim_{t\to\infty}\theta_{jk}(t)(\rho_{j}(t)-\rho_{k}(t))^{2}=0$

and this again implies either

$\lim_{t\to\infty}\min(\rho_{j}(t),\rho_{k}(t))=0\quad\text{or}\quad\lim_{t\to\infty}(\rho_{j}(t)-\rho_{k}(t))=0.$

We also use

$\sum_{i=1}^{n}\rho_{i}=1$

to get that there exists at least one $1\leq j\leq n$ such that $\rho_{j}$
does not converges to zero. The set of equilibrium of this system can be
expressed as follows:

$\mathcal{E}=\bigcup_{m=1}^{n}\mathcal{E}_{m},$

where $\mathcal{E}_{n}$ consists $(\rho_{1},\cdots,\rho_{n})$ if there exist
distinct indices $i_{1},i_{2},\cdots,i_{m}$ such that
$\rho_{i_{1}}=\cdots=\rho_{i_{m}}=\frac{1}{m}$ and other $\rho_{j}$ are zero.
From direct calculations, we have the following lemma.

######  Lemma 3.1.

Let $\rho$ be a solution to (3.5) with initial data $\rho(0)$. If
$\rho_{i}(0)=\rho_{j}(0)$, then $\rho_{i}(t)=\rho_{j}(t)$ for all $t\geq 0$.

Now, we assume the following additional condition:

$\bullet$ $\theta_{ij}$ is an increasing function of
$\min(\rho_{i},\rho_{j})$, i.e.
$\theta(\rho_{i},\rho_{j})=\tilde{\theta}(\min(\rho_{i},\rho_{j}))$ for some
$\tilde{\theta}$ and satisfies

$\min(\rho_{i},\rho_{j})\geq\min(\rho_{k},\rho_{l})\quad\Longrightarrow\quad\theta_{ij}\geq\theta_{kl}.$

As examples, we have
$\theta_{ij}=\theta(\rho_{i},\rho_{j})=\mathrm{min}(\rho_{i},\rho_{j})$. This
$\theta_{ij}$ satisfies all conditions that we assumed. From the additional
condition, we have the following property.

######  Lemma 3.2.

Suppose that the initial data $\\{\rho_{j}(0)\\}_{j=1}^{n}$ satisfy

$\rho_{1}(0)=\max_{1\leq j\leq
n}\rho_{j}(0)\quad\text{and}\quad\rho_{1}(0)>\max_{1<j\leq n}\rho_{j}(0),$

and let $\rho$ be a solution to (3.5) with initial data $\rho(0)$. Then

$\rho_{1}(t)-\max_{1<j\leq n}\rho_{j}(t)\geq\rho_{1}(0)-\max_{1<j\leq
n}\rho_{j}(0).$

###### Proof.

Without loss of generality, we set

$\rho_{2}(t)=\max_{1<j\leq
n}\rho_{j}(t)\quad\text{and}\quad\rho_{1}(t)>\rho_{2}(t)\quad\forall~{}t\in[0,\epsilon)$

for a positive $\epsilon$. From simple calculations, we have

$\displaystyle\dot{\rho}_{1}$
$\displaystyle=\kappa\sum_{k=1}^{n}\theta_{1k}(\rho_{1}-\rho_{k})=\kappa\sum_{k=1}^{n}\tilde{\theta}(\rho_{k})(\rho_{1}-\rho_{k})$
$\displaystyle\dot{\rho}_{2}$
$\displaystyle=\kappa\sum_{k=1}^{n}\theta_{2k}(\rho_{2}-\rho_{k})=\kappa\tilde{\theta}(\rho_{2})(\rho_{2}-\rho_{1})+\kappa\sum_{k=2}^{n}\tilde{\theta}(\rho_{k})(\rho_{2}-\rho_{k}).$

This implies

$\displaystyle\frac{d}{dt}(\rho_{1}-\rho_{2})$
$\displaystyle=2\kappa\tilde{\theta}(\rho_{2})(\rho_{1}-\rho_{2})+\kappa\sum_{k=2}^{n}\tilde{\theta}(\rho_{k})(\rho_{1}-\rho_{2})$
$\displaystyle=\kappa(\rho_{1}-\rho_{2})\left(2\tilde{\theta}(\rho_{2})+\sum_{k=1}^{n}\tilde{\theta}(\rho_{k})\right)\geq
0.$

So, we have

$\rho_{1}(t)-\rho_{2}(t)\geq\rho_{1}(0)-\rho_{2}(0).$

This result directly implies the desired result. ∎

Using a similar argument, we can prove the following lemma.

######  Lemma 3.3.

Suppose that the initial data $\rho(0)$ satisfy

$\rho_{1}(0)=\rho_{2}(0)=\cdots=\rho_{m}(0)>\rho_{m+1}(0)\geq\cdots\geq\rho_{n}(0),$

and let $\rho$ be a solution to (3.5) with initial data $\rho(0)$. Then, we
have

$\rho_{1}(t)=\rho_{2}(t)=\cdots=\rho_{m}(t)\quad\forall~{}t\geq 0$

and

(3.6) $\displaystyle\rho_{1}(t)-\max_{m<j\leq
n}\rho_{j}(t)\geq\rho_{1}(0)-\max_{m<j\leq n}\rho_{j}(0).$

From the previous lemma, we know that $\rho_{i}$ and $\rho_{j}$ cannot
converge to same value when $i\leq m<j$. So, we can obtain the following
theorem.

######  Theorem 3.1.

Suppose that the initial data $\rho(0)$ satisfy

$m=\left|\mathcal{S}_{M}\right|,\quad\mathcal{S}_{M}=\left\\{j:\rho_{j}(0)=\max_{1\leq
k\leq n}\rho_{k}(0),\quad 1\leq j\leq n\right\\}$

and let $\rho$ be a solution to (3.5) with initial data $\rho(0)$. Then, we
have

$\lim_{t\to\infty}\rho_{j}(t)=\frac{1}{m}\quad\forall
j\in\mathcal{S}_{M}\quad\text{and}\quad\lim_{t\to\infty}\rho_{k}(t)=0\quad\forall
j\not\in\mathcal{S}_{M}.$

###### Proof.

Without loss of generality, we set

$\rho_{1}(0)=\rho_{2}(0)=\cdots=\rho_{m}(0)>\rho_{m+1}(0)\geq\cdots\geq\rho_{n}(0).$

From the previous lemma, we know that

$\rho_{1}(t)=\rho_{2}(t)=\cdots=\rho_{m}(t)\quad\forall~{}t\geq 0$

From (3.6), we know that

$\lim_{t\to\infty}\rho_{1}(t)-\rho_{j}(t)\neq 0$

for all $j>m$. From the form of equilibrium, the limit of $\rho_{j}(t)$ should
be zero for all $j>m$. Finally, we can obtain

$\lim_{t\to\infty}\rho_{j}(t)=\frac{1}{m}$

for all $1\leq j\leq m$. ∎

######  Remark 3.1.

$\bullet$ The initial data $\rho^{0}$ is given on
$\mathcal{P}:=\\{\rho:\rho_{1}+\cdots+\rho_{n}=1,\quad\rho_{i}\geq
0\quad\forall~{}1\leq i\leq n\\}$. If we assume that the initial data is
randomly distributed on $\mathcal{P}$ uniformly then we know that the
probability of $m\geq 2$ is zero, since $\mathcal{P}$ is closed region of the
$d+1$-dimensional plane and the case $m=k$ is a part of $n-k$-dimensional
plane. Since $\mathrm{Prob}(m=1)=1$, for the generic initial data $\rho(0)$,
the complete consensus exhibits. i.e. there exists only one index $1\leq j\leq
n$ such that

$\lim_{t\to\infty}\rho_{j}(t)=1\quad\text{and}\quad\lim_{t\to\infty}\rho_{k}(t)=0\quad\forall~{}k\neq
j.$

$\bullet$ $\mathcal{E}_{1}$ is the set of stable equilibrium. Let
$\rho^{\infty}:=(\rho_{1}^{\infty},\rho_{2}^{\infty},\cdots,\rho_{n}^{\infty})=(1,0,\cdots,0)$
be an element of $\mathcal{E}_{1}$. Then, for any sufficiently small non-
negative constants $\epsilon_{2},\cdots,\epsilon_{n}$, we know that
$(\rho_{1},\rho_{2},\cdots,\rho_{n})=(1-\epsilon_{2}-\cdots-\epsilon_{n},\epsilon_{2},\cdots,\epsilon_{n})$
converges to $\rho^{\infty}$. This implies that any state
$\rho^{\infty}\in\mathcal{E}_{1}$ is a stable state.

$\bullet$ $\bigcup_{m=2}^{n}\mathcal{E}_{m}$ is the set of unstable
equilibrium. Let $2\leq k\leq n$ and
$\rho^{\infty}:=(\rho_{1}^{\infty},\rho_{2}^{\infty},\cdots,\rho_{n}^{\infty})=\Big{(}\underbrace{\frac{1}{k},\cdots,\frac{1}{k}}_{k\text{
times}},\underbrace{0,\cdots,0}_{d-k+2\text{ times}}\Big{)}$ be an element of
$\mathcal{E}_{k}$. For any constants $0<\epsilon<\frac{1}{k}$,
$(\rho_{1}^{\infty}-\epsilon,\rho_{2}^{\infty},\cdots,\rho_{d+1}^{\infty},\rho_{n}^{\infty}+\epsilon)$
does not converges to $\rho^{\infty}$. This implies that any state
$\rho^{\infty}\in\bigcup_{m=2}^{n}\mathcal{E}_{m}$ is an unstable state.

Now, we investigate the convergence rate of system when
$\theta(\rho_{i},\rho_{j})=\Big{(}\min(\rho_{i},\rho_{j})\Big{)}^{\alpha}$ for
some $\alpha\geq 1$. Obviously, this function $\theta_{ij}$ satisfies the
suggested assumptions. We assume that for a given initial configuration
$\rho^{0}$, the corresponding solution converges to
$\rho^{\infty}\in\mathcal{E}_{1}$. Without loss of generality, we assume that

$\lim_{t\to\infty}\rho_{1}(t)=1\quad\text{and}\quad\lim_{t\to\infty}\rho_{j}(t)=0\quad\forall~{}2\leq
j\leq n.$

Then we know that

$\rho_{1}(t)=\max_{1\leq i\leq
n}\rho_{i}(t),\quad\rho_{1}(t)>\frac{1}{n}\quad\forall~{}t\geq 0,$

and this yields the following calculation:

(3.7) $\displaystyle\begin{aligned}
\frac{d}{dt}\rho_{1}&=\sum_{k=1}^{n}\rho_{k}^{\alpha}(\rho_{1}-\rho_{k})\leq\sum_{k=1}^{n}(1-\rho_{1})^{\alpha}(\rho_{1}-\rho_{k})\\\
&=(1-\rho_{1})^{\alpha}((n)\rho_{1}-1)\leq(d+1)(1-\rho_{1})^{\alpha}.\end{aligned}$

Now, we have

$\displaystyle\frac{d}{dt}\rho_{1}$
$\displaystyle=\sum_{k=1}^{n}\rho_{k}^{\alpha}(\rho_{1}-\rho_{k})\geq\left(\rho_{1}-\max_{2\leq
j\leq n}\rho_{j}\right)\sum_{k=2}^{n}\rho_{k}^{\alpha}.$

From the Hölder inequality, we get

$\left(\sum_{k=2}^{n}1^{\alpha_{*}}\right)^{\frac{1}{\alpha_{*}}}\left(\sum_{k=2}^{n}\rho_{k}^{\alpha}\right)^{\frac{1}{\alpha}}\geq\left(\sum_{k=2}^{n}\rho_{k}\right),$

where $\alpha_{*}$ is the conjugate exponent of $\alpha$, i.e.
$\alpha_{*}=\frac{\alpha}{\alpha-1}$. So, we have

$\sum_{k=2}^{n}\rho_{k}^{\alpha}\geq\frac{1}{(d+1)^{\alpha-1}}\left(\sum_{k=2}^{n}\rho_{k}\right)^{\alpha}=(d+1)^{1-\alpha}(1-\rho_{1})^{\alpha}.$

Finally, we get

(3.8) $\displaystyle\begin{aligned}
\frac{d}{dt}\rho_{1}&\geq\left(\rho_{1}-\max_{2\leq j\leq
n}\rho_{j}\right)(d+1)^{1-\alpha}(1-\rho_{1})^{\alpha}\\\
&\geq\left(\rho_{1}^{0}-\max_{2\leq j\leq
n}\rho_{j}^{0}\right)(d+1)^{1-\alpha}(1-\rho_{1})^{\alpha}.\end{aligned}$

Here, we use Lemma 3.2 in the last inequality. Now, we combine (3.7) and (3.8)
to get

(3.9) $\displaystyle C_{1}(1-\rho_{1})^{\alpha}\leq\frac{d}{dt}\rho_{1}\leq
C_{2}(1-\rho_{1})^{\alpha}$

for some positive constants $C_{1}$ and $C_{2}$. From solving the ODE, we have
the following lemma.

######  Lemma 3.4.

Let $x(t)$ be a solution to the following ODE:

$\dot{x}(t)=C(1-x(t))^{\alpha},\quad x(0)=x^{0}.$

for some positive constant $C$ and $\alpha\geq 1$. Then, we have

$x(t)=\begin{cases}1-(1-x^{0})e^{-Ct}&\quad\text{when }\alpha=1,\\\
\displaystyle
1-\frac{1}{\left((1-x^{0})^{1-\alpha}+C(\alpha-1)t\right)^{\frac{1}{\alpha-1}}}&\quad\text{when
}\alpha>1.\end{cases}$

We combine (3.9) and Lemma 3.4 to get

$(1-\rho_{1}^{0})e^{-C_{1}t}\leq
1-\rho_{1}\leq(1-\rho_{1}^{0})e^{-C_{2}t}\quad\text{when }\alpha=1$

and

$\frac{1}{\left((1-\rho_{1}^{0})^{1-\alpha}+C_{1}(\alpha-1)t\right)^{\alpha-1}}\leq
1-\rho_{1}\leq\frac{1}{\left((1-\rho_{1}^{0})^{1-\alpha}+C_{2}(\alpha-1)t\right)^{\alpha-1}}\quad\text{when
}\alpha>1.$

So, we can conclude the result as follows:

(3.10) $\displaystyle 1-\rho_{1}(t)\simeq\begin{cases}\text{blow up in a
finite time}\quad&\text{when}\quad 0<\alpha<1,\\\ \text{converges to zero
exponentially}&\text{when}\quad\alpha=1,\\\
Ct^{(1-\alpha)^{-1}}&\text{when}\quad 1<\alpha,\\\ \end{cases}$

For example, when $\alpha=1,2,3$, we have

(3.11) $\displaystyle 1-\rho_{1}(t)\simeq e^{-Ct},\quad
1-\rho_{1}(t)\simeq\frac{1}{t},\quad 1-\rho_{1}(t)\simeq\frac{1}{\sqrt{t}},$

respectively. So we can conclude that the convergence rate is sensitive to
$\alpha\geq 1$. So far, we investigated emergent behaviors on the complete
graph. From now on, we consider the system (3.1) on general graphs. We set the
network topology as follows:

(3.12) $\displaystyle\omega_{ij}=\begin{cases}1\quad\text{when }(i,j)\in E,\\\
0\quad\text{when }(i,j)\not\in E.\end{cases}$

From a simple calculation, we can obtain the following system:

(3.13) $\displaystyle\frac{d\rho_{j}}{dt}=\kappa\sum_{l\in
N_{j}}\theta_{jl}(\rho_{j}-\rho_{l}).$

This yields

$\displaystyle\frac{d}{dt}\sum_{j=1}^{n}\rho_{j}^{2}$
$\displaystyle=\sum_{j=1}^{n}2\kappa\rho_{j}\left(\sum_{l\in
N_{j}}\theta_{jl}(\rho_{j}-\rho_{l})\right)$
$\displaystyle=2\kappa\sum_{(j,l)\in
E}\theta_{jl}\rho_{j}(\rho_{j}-\rho_{l})=\kappa\sum_{(j,l)\in
E}\theta_{jl}(\rho_{j}-\rho_{l})^{2}\geq 0.$

From the continuity of $\theta_{jl}$, we also apply Barbalat’s lemma [1] to
get

$\displaystyle\lim_{t\to\infty}\sum_{(j,l)\in
E}\theta_{jl}(\rho_{j}-\rho_{l})^{2}=0.$

Since each term in the sum is non-negative, we again obtain

$\lim_{t\to\infty}(\rho_{j}-\rho_{l})^{2}\theta_{jl}=0\quad\forall~{}(j,l)\in
E.$

From this relation, if $(j,l)\in E$ then we have either

(3.14)
$\displaystyle\lim_{t\to\infty}(\rho_{j}-\rho_{l})=0\quad\text{or}\quad\lim_{t\to\infty}\min(\rho_{j},\rho_{l})=0.$

When the graph is complete(i.e. $E=(V\times V)\backslash\mathrm{diag}(V\times
V)$), then (3.14) holds for all $j,l\in V$ with $j\neq l$. This yields that if
there exists one index $i\in V$ such that $\rho_{i}(0)>\max_{j\neq
i}\rho_{j}(0)$, then

(3.15) $\displaystyle\rho_{i}\to 1,\quad\rho_{j}\to 0\quad\forall~{}j\in
V\backslash\\{i\\}.$

We have already shown the above property in Theorem 3.1. However, when the
graph is not a complete graph, then we have the counter example for (3.15). We
consider the graph is given as a square. We denote four vertices as $A$, $B$,
$C$ and $D$, counter-clockwise order.

We also assume that the initial density is given as

$\rho_{A}(0)=0.6,\quad\rho_{B}(0)=0.1,\quad\rho_{C}(0)=0.2,\quad\rho_{D}(0)=0.1.$

Then, the we can easily check that

$\rho_{A}(t)\text{ and }\rho_{C}(t)\quad\text{increase},\quad\rho_{B}(t)\text{
and }\rho_{D}(t)\quad\text{decrease}.$

Using (3.15), we can conclude that

$\lim_{t\to\infty}\rho_{A}(t)=\rho_{A}^{\infty},\quad\lim_{t\to\infty}\rho_{C}(t)=\rho_{C}^{\infty}$

for some $\rho_{A}^{\infty}>0.6$ and $\rho_{C}^{\infty}>0.2$ which satisfies
$\rho_{A}^{\infty}+\rho_{C}^{\infty}=1$, and

$\lim_{t\to\infty}\rho_{B}(t)=\lim_{t\to\infty}\rho_{D}(t)=0.$

We can say that the final state $(\rho_{A}^{\infty},0,\rho_{C}^{\infty},0)$
corresponds to the bipolar state of the Kuramoto model. We can summarize the
above results as the following proposition.

######  Proposition 3.1.

Let $(\rho,S)$ be a solution to system (3.19), and the network topology is
given as (3.12). If $(j,l)\in E$, then (3.14) holds.

###  3.3. Second order dynamics on the two points graph

In this subsection, we study emergent behavior of system (3.3) with general
initial configuration when $n=2$. System (3.3) with $n=2$ can be written as

$\left\\{\begin{aligned} \frac{d\rho_{1}}{dt}&=(S_{1}-S_{2})\theta_{12},\\\
\frac{d\rho_{2}}{dt}&=-(S_{1}-S_{2})\theta_{12},\\\
\frac{dS_{1}}{dt}&=\kappa^{2}(\rho_{1}-\rho_{2})\theta_{12}+\frac{1}{2}(-(S_{1}-S_{2})^{2}+\kappa^{2}(\rho_{1}-\rho_{2})^{2})\frac{\partial\theta_{12}}{\partial\rho_{1}},\\\
\frac{dS_{2}}{dt}&=-\kappa^{2}(\rho_{1}-\rho_{2})\theta_{12}+\frac{1}{2}(-(S_{1}-S_{2})^{2}+\kappa^{2}(\rho_{1}-\rho_{2})^{2})\frac{\partial\theta_{12}}{\partial\rho_{2}}.\end{aligned}\right.$

This yields

$\left\\{\begin{aligned}
&\frac{d}{dt}(\rho_{1}-\rho_{2})=2(S_{1}-S_{2})\theta_{12},\\\
&\frac{d}{dt}(S_{1}-S_{2})=2\kappa^{2}(\rho_{1}-\rho_{2})\theta_{12}-((S_{1}-S_{2})^{2}-\kappa^{2}(\rho_{1}-\rho_{2})^{2})\left(\frac{\partial\theta_{12}}{\partial\rho_{1}}-\frac{\partial\theta_{12}}{\partial\rho_{2}}\right).\end{aligned}\right.$

Now, we substitute $\rho_{1}=r$, $\rho_{2}=1-r$, $S=S_{1}-S_{2}$, and
$\theta_{12}=\theta$ into the above system to get

(3.16) $\left\\{\begin{aligned} &\frac{dr}{dt}=S\theta,\\\
&\frac{dS}{dt}=2\kappa^{2}(2r-1)\theta-\left(S^{2}-\kappa^{2}(2r-1)^{2}\right)\frac{d\theta}{dr}.\end{aligned}\right.$

Here we used

$\frac{d\theta}{dr}=\frac{\partial\theta_{12}}{\partial\rho_{1}}\frac{d\rho_{1}}{dr}+\frac{\partial\theta_{12}}{\partial\rho_{2}}\frac{d\rho_{2}}{dr}=\frac{\partial\theta_{12}}{\partial\rho_{1}}-\frac{\partial\theta_{12}}{\partial\rho_{2}}.$

Recall that the following Hamiltonian

(3.17) $\displaystyle\mathcal{H}=\frac{\theta}{2}(S^{2}-\kappa^{2}(2r-1)^{2})$

is a constant of motion. If $\mathcal{H}(0)=:\mathcal{H}_{0}=0$, then it can
be reduced to a gradient flow formulation(subsection 3.2). Let assume that
$\mathcal{H}_{0}>0$. Then, we get

$S^{2}=\kappa^{2}(2r-1)^{2}+\frac{2\mathcal{H}_{0}}{\theta}>0.$

This implies that $S(t)\neq 0$ for all $t$. From this fact, if we assume the
initial data satisfies $S_{0}>0$ then $S(t)>0$ for all $t$. This yields

$S=\sqrt{\kappa^{2}(2r-1)^{2}+\frac{2\mathcal{H}_{0}}{\theta}}.$

Since $S\theta>0$ for all $t$, we know that $r$ increases. Also, we know that
$S\geq\sqrt{\frac{2\mathcal{H}_{0}}{\theta_{0}}}>0$. This implies that

$\lim_{t\to\infty}r(t)=1.$

Now, we are interested in the convergence rate of $r$. We can easily deduce
that $S$ tends to infinity as $r\to 1$, since $\theta\big{|}_{r=1}=0$. Since
$\kappa^{2}(2r-1)^{2}$ is bounded, we can say that

$S\simeq\sqrt{\frac{2\mathcal{H}_{0}}{\theta}}.$

Then, we have

$\frac{dr}{dt}=S\theta\simeq\sqrt{2\mathcal{H}_{0}\theta}.$

Since $r$ converges to one, we assume have that $\min(r,1-r)=1-r$. We assume
that $\theta\simeq(1-r)^{\alpha}$ for some $\alpha$. This yields

$\frac{\dot{r}}{(1-r)^{\alpha/2}}\simeq\sqrt{2\mathcal{H}_{0}}$

If $\alpha=2$, we get

$\ln(1-r_{0})-\ln(1-r(t))\simeq\sqrt{2\mathcal{H}_{0}}t,$

and this implies

$1-r(t)\simeq(1-r_{0})e^{-\sqrt{2\mathcal{H}_{0}}t}.$

If $\alpha>2$, then we have

$\left[-\left(1-r\right)^{1-\alpha/2}\right]_{0}^{t}\simeq\left(\frac{\alpha}{2}-1\right)\sqrt{2\mathcal{H}_{0}}t.$

This yields

$(1-r(t))^{1-\alpha/2}\simeq(1-r_{0})^{1-\alpha/2}+\left(\frac{\alpha}{2}-1\right)\sqrt{2\mathcal{H}_{0}}t,$

or equivalently

$1-r(t)\simeq\frac{1}{\left((1-r_{0})^{1-\alpha/2}+\left(\frac{\alpha}{2}-1\right)\sqrt{2\mathcal{H}_{0}}t\right)^{(\alpha/2-1)^{-1}}}\simeq\frac{C}{t^{(\alpha/2-1)^{-1}}}.$

If $0<\alpha<2$, then the solution blow up in a finite time. We can summarize
the above results as follows:

(3.18) $\displaystyle 1-r(t)\simeq\begin{cases}\text{blow up in a finite
time}\quad&\text{when}\quad 0<\alpha<2,\\\
Ce^{-\sqrt{2\mathcal{H}_{0}}t}&\text{when}\quad\alpha=2,\\\
Ct^{(1-\alpha/2)^{-1}}\quad&\text{when}\quad
2<\alpha,\end{cases}\quad\text{as}\quad t\to\infty.$

We can summarize the above result below.

######  Proposition 3.2.

Let $(r,S)$ be a solution of (3.16) with the initial data $(r_{0},S_{0})$
which satisfies the following:

$\mathcal{H}_{0}=\mathcal{H}(0)>0,\quad S_{0}>0,$

where $\mathcal{H}$ is defined in (3.17). Then, $r(t)$ satisfies (3.18).

######  Remark 3.2.

We can compare two results (3.10) and (3.18). The case (3.10) occurs when
$\mathcal{H}_{0}=0$ and the case (3.18) occurs when $\mathcal{H}_{0}>0$. So,
we can check the existence of bifurcation at $\mathcal{H}_{0}=0$.

We will check general case $n>2$ in the numerical example section (Section
5.3).

###  3.4. Hopf-Cole transform on discrete graphs

In this subsection, we reduce the following system (2.7) onto discrete graphs
and provide its Hopf-Cole transform. For the notation simplicities, we set

$\partial_{k}\mathcal{F}=\frac{\partial\mathcal{F}}{\partial\rho_{k}},\quad\partial^{2}_{kl}\mathcal{F}=\frac{\partial^{2}\mathcal{F}}{\partial\rho_{k}\partial\rho_{l}}\quad\forall~{}1\leq
k,l\leq n.$

From a similar argument which used in [5], we can discretize system (2.7) as
follows:

(3.19) $\displaystyle\begin{cases}\displaystyle\frac{d\rho_{j}}{dt}-\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-S_{l})\theta_{jl}=0,\vspace{0.2cm}\\\
\displaystyle\frac{dS_{j}}{dt}+\frac{1}{2}\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-S_{l})^{2}\frac{\partial\theta_{jl}}{\partial\rho_{j}}=\frac{\partial}{\partial\rho_{j}}\left(\frac{1}{2}\sum_{(k,l)\in
E}\frac{\omega_{kl}}{2}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})^{2}\theta_{kl}\right).\end{cases}$

Also, the initial condition of continuous domain version system (2.7)

$S_{0}=-\frac{\delta\mathcal{F}(\rho_{0})}{\delta\rho_{0}}$

can be reduced as follows:

(3.20) $\displaystyle
S_{j}(0)=-\partial_{j}\mathcal{F}(0):=-\frac{\partial\mathcal{F}(\rho)}{\partial\rho_{j}}\Bigg{|}_{\rho=\rho(0)},\quad\rho=(\rho_{1},\cdots,\rho_{n}).$

Here, $\mathcal{F}$ is a kind of potential functional of the system. When the
initial condition of system (3.19) is given as (3.20), then the system can be
reduce to the gradient flow system with potential $\mathcal{F}$. On the other
hand, if the initial condition satisfies

(3.21) $\displaystyle
S_{j}(0)=-\partial_{j}\mathcal{F}(0):=-\frac{\partial\mathcal{F}(\rho)}{\partial\rho_{j}}\Bigg{|}_{\rho=\rho(0)},\quad\rho=(\rho_{1},\cdots,\rho_{n}).$

then the system can be reduced to the gradient flow system with potential
$-\mathcal{F}$. This phenomena can be interpreted as that initial conditions
(3.20) and (3.21) yield a gradient descent flow and a gradient upscent flow,
respectively.

We can also express system (3.19) as follows:

(3.22)
$\displaystyle\begin{cases}\displaystyle\frac{d}{dt}\rho_{j}=\frac{\partial}{\partial
S_{j}}\mathcal{H}(\rho,S),\vspace{0.2cm}\\\
\displaystyle\frac{d}{dt}S_{j}=-\frac{\partial}{\partial\rho_{j}}\mathcal{H}(\rho,S),\end{cases}$

where the Hamiltonian is given as follows:

$\mathcal{H}(\rho,S)=\frac{1}{4}\sum_{(i,j)\in
E}\omega_{ij}\theta_{ij}\Big{(}(S_{i}-S_{j})^{2}-(\partial_{i}\mathcal{F}-\partial_{j}\mathcal{F})^{2}\Big{)}.$

So we can conclude that the discretized system (3.19) preserve the property of
Hamiltonian systems, and this means the system is well-discretized. From
relation (3.22), we can easily check that

$\frac{d}{dt}\mathcal{H}(\rho(t),S(t))=0$

when $(\rho(t),S(t))$ is a solution to system (3.19).

Now, we want to obtain the discretized version of the generalized Hopf-Cole
transformation (2.10). To mimic substitution (2.9) of the continuous version,
we consider the the following substitution:

(3.23)
$\displaystyle\begin{cases}\partial_{j}\mathcal{F}=\xi_{j}+\xi^{*}_{j},\\\
S_{j}=\xi_{j}-\xi^{*}_{j},\end{cases}\quad\Longleftrightarrow\qquad\begin{cases}\xi_{j}=\frac{1}{2}(\partial_{j}\mathcal{F}+S_{j}),\\\
\xi^{*}_{j}=\frac{1}{2}(\partial_{j}\mathcal{F}-S_{j}).\end{cases}$

Then, we have the following proposition.

######  Proposition 3.3 (Hopf-Cole transform on a finite graph).

Let $(\rho,S)$ be a solution to (3.19), and the substitution $(\xi,\xi^{*})$
is given as (3.23). Then $(\xi,\xi^{*})$ follows the following dynamics:

(3.24)
$\displaystyle\begin{cases}\displaystyle\frac{d\xi_{j}}{dt}=\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(\xi_{k}-\xi_{l})+\sum_{l\in
N_{j}}\omega_{jl}(\xi_{l}^{*}-\xi_{j}^{*})(\xi_{l}-\xi_{j})\frac{\partial\theta_{jl}}{\partial\rho_{j}},\\\
\displaystyle\frac{d\xi^{*}_{j}}{dt}=-\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(\xi_{k}^{*}-\xi_{l}^{*})-\sum_{l\in
N_{j}}\omega_{jl}(\xi_{l}^{*}-\xi_{j}^{*})(\xi_{l}-\xi_{j})\frac{\partial\theta_{jl}}{\partial\rho_{j}}.\end{cases}$

###### Proof.

From (3.19) and (3.23), we have the following calculation:

$\displaystyle\frac{d\xi_{j}}{dt}$
$\displaystyle=\frac{1}{2}\frac{d\partial_{j}\mathcal{F}}{dt}+\frac{1}{2}\frac{dS_{j}}{dt}$
$\displaystyle=\frac{1}{2}\sum_{k}\left(\partial^{2}_{jk}\mathcal{F}\frac{d\rho_{k}}{dt}\right)+\frac{1}{2}\frac{dS_{j}}{dt}$
$\displaystyle=\frac{1}{2}\sum_{\begin{subarray}{c}k,\\\ l\in
N_{k}\end{subarray}}\partial^{2}_{jk}\mathcal{F}\omega_{kl}(S_{k}-S_{l})\theta_{kl}-\frac{1}{4}\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-S_{l})^{2}\frac{\partial\theta_{jl}}{\partial\rho_{j}}$
$\displaystyle\quad+\frac{1}{4}\sum_{(k,l)\in
E}\frac{\omega_{kl}}{2}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})^{2}\frac{\partial\theta_{kl}}{\partial\rho_{j}}+\frac{1}{4}\sum_{(k,l)\in
E}\omega_{kl}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})(\partial^{2}_{kj}\mathcal{F}-\partial^{2}_{lj}\mathcal{F})\theta_{kl}$
$\displaystyle=\frac{1}{2}\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(S_{k}-S_{l})+\frac{1}{2}\sum_{(k,l)\in
E}\partial^{2}_{kj}\mathcal{F}\omega_{kl}\theta_{kl}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})$
$\displaystyle\quad-\frac{1}{4}\sum_{l\in
N_{j}}\omega_{jl}\Big{(}(S_{j}-S_{l})^{2}-(\partial_{k}\mathcal{F}-\delta_{l}\mathcal{F})^{2}\Big{)}\frac{\partial\theta_{jl}}{\partial\rho_{j}}$
$\displaystyle=-\frac{1}{2}\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(-S_{k}-\partial_{k}\mathcal{F}+S_{l}+\partial_{l}\mathcal{F})$
$\displaystyle\quad-\frac{1}{4}\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-\partial_{j}\mathcal{F}-S_{l}+\partial_{l}\mathcal{F})(S_{j}+\partial_{j}\mathcal{F}-S_{l}-\partial_{l}\mathcal{F})\frac{\partial\theta_{jl}}{\partial\rho_{j}}$
$\displaystyle=\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(\xi_{k}-\xi_{l})+\sum_{l\in
N_{j}}\omega_{jl}(\xi_{l}^{*}-\xi_{j}^{*})(\xi_{l}-\xi_{j})\frac{\partial\theta_{jl}}{\partial\rho_{j}}.$

From a similar calculation, we have

$\displaystyle\frac{d\xi^{*}_{j}}{dt}$
$\displaystyle=\frac{1}{2}\frac{d\partial\mathcal{F}_{j}}{dt}-\frac{1}{2}\frac{dS_{j}}{dt}$
$\displaystyle=\frac{1}{2}\sum_{k}\left(\partial^{2}_{jk}\mathcal{F}\frac{d\rho_{k}}{dt}\right)-\frac{1}{2}\frac{dS_{j}}{dt}$
$\displaystyle=\frac{1}{2}\sum_{\begin{subarray}{c}k,\\\ l\in
N_{k}\end{subarray}}\partial^{2}_{jk}\mathcal{F}\omega_{kl}(S_{k}-S_{l})\theta_{kl}+\frac{1}{4}\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-S_{l})^{2}\frac{\partial\theta_{jl}}{\partial\rho_{j}}$
$\displaystyle\quad-\frac{1}{4}\sum_{(k,l)\in
E}\frac{\omega_{kl}}{2}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})^{2}\frac{\partial\theta_{kl}}{\partial\rho_{j}}-\frac{1}{4}\sum_{(k,l)\in
E}\omega_{kl}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})(\partial^{2}_{kj}\mathcal{F}-\partial^{2}_{lj}\mathcal{F})\theta_{kl}$
$\displaystyle=\frac{1}{2}\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(S_{k}-S_{l})-\frac{1}{2}\sum_{(k,l)\in
E}\partial^{2}_{kj}\mathcal{F}\omega_{kl}\theta_{kl}(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})$
$\displaystyle\quad+\frac{1}{4}\sum_{l\in
N_{j}}\omega_{jl}\Big{(}(S_{j}-S_{l})^{2}-(\partial_{k}\mathcal{F}-\partial_{l}\mathcal{F})^{2}\Big{)}\frac{\partial\theta_{jl}}{\partial\rho_{j}}$
$\displaystyle=-\frac{1}{2}\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(-S_{k}+\partial_{k}\mathcal{F}+S_{l}-\partial_{l}\mathcal{F})$
$\displaystyle\quad+\frac{1}{4}\sum_{l\in
N_{j}}\omega_{jl}(S_{j}-\partial_{j}\mathcal{F}-S_{l}+\partial_{l}\mathcal{F})(S_{j}+\partial_{j}\mathcal{F}-S_{l}-\partial_{l}\mathcal{F})\frac{\partial\theta_{jl}}{\partial\rho_{j}}$
$\displaystyle=-\sum_{(k,l)\in
E}\partial^{2}_{jk}\mathcal{F}\omega_{kl}\theta_{kl}(\xi_{k}^{*}-\xi_{l}^{*})-\sum_{l\in
N_{j}}\omega_{jl}(\xi_{l}^{*}-\xi_{j}^{*})(\xi_{l}-\xi_{j})\frac{\partial\theta_{jl}}{\partial\rho_{j}}.$

So, we have the desired result. ∎

If the initial data $(\rho(0),S(0))$ of system (3.19) satisfies (3.20), then
the initial data $(\xi,\xi^{*})$ of system (3.24) should satisfy

(3.25) $\displaystyle\xi_{j}(0)=0\quad\forall~{}j=1,2,\cdots,n.$

Here, $n$ is the number of vertices.

######  Lemma 3.5.

Let $(\xi,\xi^{*})$ be a solution to system (3.24) with the initial data
$(\xi(0),\xi^{*}(0))$. If the initial data satisfies (3.25), we have

$\xi_{j}(0)=0\quad\forall~{}j=1,2,\cdots,n,\quad t\geq 0.$

Furthermore, $\xi^{*}$ follows the following dynamics:

$\frac{d}{dt}\xi_{j}^{*}=-\sum_{(k,l)\in
E}\partial^{2}\mathcal{F}_{jk}\theta_{kl}(\xi_{k}^{*}-\xi_{l}^{*}).$

###### Proof.

From the first equation of system (3.24), we can express

$\frac{d\xi}{dt}(t)=A(t)\xi(t)$

for some time-dependent matrix of size $n\times n$. From the uniqueness of
system, we can obtain the desired result. ∎

######  Remark 3.3.

If we discretize system (3.19), then such numeric scheme may not preserve
$S_{j}+\delta\mathcal{F}_{j}$ due to the numeric errors. However, if we use
(3.24) to discretize system, then $\xi_{j}(0)\equiv 0$ for all $j$ guarantee
that $\xi_{j}\equiv 0$ for all time.

Corresponding lemma for Lemma 3.5 on continuous domain can be organized as
Lemma A.3.

##  4\. Two points graph: analytical solutions.

In this section, we study system (3.19) with $n=2$ analytically.

System (3.19) with $n=2$ can be written as

(4.1)
$\displaystyle\begin{cases}\displaystyle\frac{d\rho_{1}}{dt}-(S_{1}-S_{2})\theta_{12}=0,\vspace{0.2cm}\\\
\displaystyle\frac{d\rho_{2}}{dt}-(S_{2}-S_{1})\theta_{12}=0,\vspace{0.2cm}\\\
\displaystyle\frac{dS_{1}}{dt}+\frac{1}{2}(S_{1}-S_{2})^{2}\frac{\partial\theta_{12}}{\partial\rho_{1}}=\frac{\partial}{\partial\rho_{1}}\left(\frac{1}{2}(\partial_{1}\mathcal{F}-\partial_{2}\mathcal{F})^{2}\theta_{12}\right),\vspace{0.2cm}\\\
\displaystyle\frac{dS_{2}}{dt}+\frac{1}{2}(S_{1}-S_{2})^{2}\frac{\partial\theta_{12}}{\partial\rho_{2}}=\frac{\partial}{\partial\rho_{2}}\left(\frac{1}{2}(\partial_{1}\mathcal{F}-\partial_{2}\mathcal{F})^{2}\theta_{12}\right).\end{cases}$

Recall that the Hamiltonian of system (4.1) can be expressed as

$\mathcal{H}(\rho_{1},\rho_{2},S_{1},S_{2})=\frac{1}{2}\theta_{12}((S_{1}-S_{2})^{2}-(\partial_{1}\mathcal{F}-\partial_{2}\mathcal{F})^{2}).$

If we substitute

$\rho_{1}=r,\quad\rho_{2}=1-r,\quad S-1-S_{2}=S,\quad\theta_{12}=\theta,$

then system (4.1) yields the following dynamics of $r$ and $S$:

$\begin{cases}\displaystyle\frac{dr}{dt}-S\theta=0,\vspace{0.2cm}\\\
\displaystyle\frac{dS}{dt}+\frac{1}{2}S^{2}\frac{d\theta}{dr}=\frac{d}{dr}\left(\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right)\end{cases}$

and the corresponding Hamiltonian is

$\mathcal{H}(r,S)=\frac{\theta}{2}\left(S^{2}-\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right).$

The corresponding Lagrangian is

$\mathcal{L}(r,S)=\frac{\theta}{2}\left(S^{2}+\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right).$

We first demonstrate that the Hamiltonian flow (4.1) has a variational
formulation.

######  Proposition 4.1.

Consider the following Lagrangian action minimization problem:

(4.2) $\displaystyle\mathcal{A}(r_{0},r_{1}):=\text{minimize}$
$\displaystyle\int_{0}^{1}\left(\frac{\dot{r}^{2}}{2\theta}+\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right)dt$
subject to $\displaystyle r(0)=r_{0},\quad r(1)=r_{1}.$

The critical point system of minimization problem (4.2) satisfies the
Hamiltonian flow (4.1).

###### Proof.

Since $S=\frac{\dot{r}}{\theta}$, we can rewrite the Hamiltonian as follows:

$\mathcal{H}(r,\dot{r})=\frac{\dot{r}^{2}}{2\theta}-\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2},$

where we consider the above Hamiltonian as a function of $r$ and $\dot{r}$,
since $\theta$ is a function of $r$. Since the Hamiltonian is constant along
time-evolution, we assume that

$\frac{\dot{r}^{2}}{2\theta}-\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2}=\mathcal{H}_{0}.$

For the notation simplicity, we set

$r(0)=r_{0}\quad\text{and}\quad r(1)=r_{1}.$

We introduce the following substitution:

$\frac{dx}{dr}=\frac{1}{\sqrt{\theta}}.$

Then, we get

(4.3) $\displaystyle\frac{d^{2}}{dt^{2}}x$
$\displaystyle=\frac{d}{dt}\left(\frac{dx}{dr}\dot{r}\right)=\frac{d}{dt}\frac{\dot{r}}{\sqrt{\theta}}=\frac{d}{dr}\left(\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right)\sqrt{\theta}=\frac{d}{dx}\left(\frac{1}{2}\left(\frac{d\mathcal{F}}{dx}\right)^{2}\right)$

From the definition of the action, we have

(4.4) $\displaystyle\mathcal{A}$
$\displaystyle=\int_{0}^{1}\mathcal{L}(r,\dot{r})dt=\int_{0}^{1}\left(\frac{\dot{r}^{2}}{2\theta}+\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right)dt=\int_{0}^{1}\left(\frac{1}{2}\dot{x}^{2}+\frac{1}{2}\left(\frac{d\mathcal{F}}{dx}\right)^{2}\right)dt.$

∎

By the following proposition, we can obtain the explicit form of the solution.

######  Proposition 4.2.

If $\mathcal{F}$ satisfies

(4.5)
$\displaystyle\mathcal{F}(r)=\frac{1}{2}\left(\int_{1/2}^{r}\frac{1}{\sqrt{\theta(r_{*})}}dr_{*}\right)^{2},$

then we have the following explicit form of the solution $r(t)$:

(4.6)
$\displaystyle\int_{1/2}^{r(t)}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}=\frac{1}{\sinh
1}\left(\sinh(1-t)\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}+\sinh
t\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right).$

In addition, the least action minimizer satisfies

$\mathcal{A}(r_{0},r_{1})=\frac{\cosh 1}{\sinh
1}\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}-\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\frac{2\cosh
1-1}{\sinh
1}\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}.$

###### Proof.

From the definition of $x(r)$, we get

$\mathcal{F}(r)=\frac{1}{2}x(r)^{2}.$

Since $x$ solves (4.3) and $\mathcal{F}(r)=\frac{1}{2}x(r)^{2}$, we have the
general form of $x(t)$ as follows:

$x(t)=a\cosh t+b\sinh t,$

where $a=x(0)$ and $b=\dot{x}(0)$. Since $x(0)$ and $x(1)$ are fixed, we have

$a=x(0),\quad b=\frac{x(1)-x(0)\cosh 1}{\sinh 1}.$

Then, we have

$x(t)=x(0)\cosh t+\left(\frac{x(1)-x(0)\cosh 1}{\sinh 1}\right)\sinh
t=x(0)\frac{\sinh(1-t)}{\sinh 1}+x(1)\frac{\sinh t}{\sinh 1}.$

Finally, we substitute

$x(t)=\int_{1/2}^{r(t)}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}$

into the above equation to get the desired result. From (4.6), we know that

$\frac{\dot{r}}{\sqrt{\theta(r_{*})}}=\frac{1}{\sinh
1}\left(-\cosh(1-t)\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}+\cosh
t\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)$

and

$\frac{d\mathcal{F}}{dr}=\frac{1}{\sqrt{\theta(r)}}\int^{r}_{1/2}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}=\frac{1}{\sinh
1\sqrt{\theta(r)}}\left(\sinh(1-t)\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}+\sinh
t\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right).$

We substitute (4.6) into $\mathcal{A}$ as follows:

$\displaystyle\int_{0}^{1}\left(\frac{\dot{r}^{2}}{2\theta}+\frac{\theta}{2}\left(\frac{d\mathcal{F}}{dr}\right)^{2}\right)dt$
$\displaystyle=\frac{1}{2\sinh^{2}1}\int_{0}^{1}\left(-\cosh(1-t)\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}+\cosh
t\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}dt$
$\displaystyle+\frac{1}{2\sinh^{2}1}\int_{0}^{1}\left(\sinh(1-t)\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}+\sinh
t\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}dt$
$\displaystyle=\frac{1}{2\sinh^{2}1}\left(\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\left(\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}\right)\int_{0}^{1}(\cosh^{2}(1-t)+\sinh^{2}(1-t))dt$
$\displaystyle-\frac{1}{\sinh^{2}1}\int_{0}^{1}(\cosh(1-t)\cosh
t-\sinh(1-t)\sinh
t)\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}$
$\displaystyle=\frac{\sinh
2}{2\sinh^{2}1}\left(\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\left(\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}\right)-\frac{\sinh
1}{\sinh^{2}1}\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}$
$\displaystyle=\frac{\cosh 1}{\sinh
1}\left(\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\left(\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}\right)-\frac{1}{\sinh
1}\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}$
$\displaystyle=\frac{\cosh 1}{\sinh
1}\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}-\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\frac{2\cosh
1-1}{\sinh
1}\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}.$

So, we can conclude that the

$\mathcal{A}(r_{0},r_{1})=\frac{\cosh 1}{\sinh
1}\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}-\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\frac{2\cosh
1-1}{\sinh
1}\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}.$

We can easily check that $\mathcal{A}(r_{0},r_{1})\geq 0$ for any $r_{0}$ and
$r_{1}$ and equality only holds for $r_{0}=r_{1}=\frac{1}{2}$, since

$\cosh 1x^{2}+(2\cosh 1-1)xy+\cosh 1y^{2}\geq 0$

and the equality only holds for $x=y=0$($\because D=(2\cosh
1-1)^{2}-4\cosh^{2}1<0$). ∎

######  Proposition 4.3.

We formulate a divergence function on a two point graph:

$\begin{split}D(r_{0},r_{1})&=\mathcal{A}(r_{0},r_{1})-\frac{1}{2}\mathcal{A}(r_{0},r_{0})-\frac{1}{2}\mathcal{A}(r_{1},r_{1})\\\
&=\frac{1}{2\sinh
1}\left(\int_{r_{0}}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}.\end{split}$

######  Remark 4.1.

We remark that the above divergence function $D(r_{0},r_{1})$ is a
generalization of the discrete Wasserstein metric; see related studies in AI
inference problems [4]. In our selection of $\mathcal{F}$ in (4.5), $D$
recovers a two point discrete Wasserstein distances derived in [17].

###### Proof.

The proof follows from a direct calculation.

$\displaystyle D(r_{0},r_{1})$
$\displaystyle=\mathcal{A}(r_{0},r_{1})-\frac{1}{2}\mathcal{A}(r_{0},r_{0})-\frac{1}{2}\mathcal{A}(r_{1},r_{1})$
$\displaystyle=\frac{\cosh 1}{\sinh
1}\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}-\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\frac{2\cosh
1-1}{\sinh
1}\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}$
$\displaystyle\quad-\frac{2\cosh 1-1}{2\sinh
1}\left[\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}+\left(\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}\right]$
$\displaystyle=\left(\frac{\cosh 1}{\sinh 1}-\frac{2\cosh 1-1}{2\sinh
1}\right)\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}-\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}$
$\displaystyle=\frac{1}{2\sinh
1}\left(\int_{1/2}^{r_{0}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}-\int_{1/2}^{r_{1}}\frac{dr_{*}}{\sqrt{\theta(r_{*})}}\right)^{2}.$

∎

Now, we introduce some examples $\mathcal{F}$ which exhibits aggregation
phenomena. Recall that the Kuramoto type potential (3.2) is minimized when

$\rho_{i}=\begin{cases}1\quad\text{when }i=j,\\\ 0\quad\text{when }i\neq
j\end{cases}$

for some fixed index $1\leq j\leq n$. Entropy functionals satisfy similar
property. So, we review various types of entropies.

######  Example 4.1.

Let $\rho=(\rho_{1},\cdots,\rho_{n})$ be a probability distribution given on
$n$-point graph. Then, we have the following entropies.

(1) Shannon entropy:

$\displaystyle H(\rho):=-\sum_{i=1}^{n}\rho_{i}\log\rho_{i}.$

(2) Rényi entropy:

$\displaystyle
H_{\alpha}(\rho)=\frac{1}{1-\alpha}\log\left(\sum_{i=1}^{n}\rho_{i}^{\alpha}\right),$

where $\alpha$ is a non-negative constant with $\alpha\neq 1$. From the
L’H$\hat{o}$spital theorem, we can easily check that Rényi entropy converges
to Shannon entropy as $\alpha$ goes to 1. i.e.

$\lim_{\alpha\to 1}H_{\alpha}(\rho)=H(\rho).$

Similarly, we can also check that

$\displaystyle\lim_{\alpha\to\infty}H_{\alpha}(\rho)$
$\displaystyle=\lim_{\alpha\to\infty}\frac{1}{1-\alpha}\log\left(\sum_{i=1}^{n}\rho_{i}^{\alpha}\right)=\lim_{\alpha\to\infty}\frac{\sum_{i=1}^{n}\log(\rho_{i})\rho_{i}^{\alpha}}{-\sum_{i=1}^{n}\rho_{i}^{\alpha}}$
$\displaystyle=-\lim_{\alpha\to\infty}\frac{\sum_{i=1}^{n}\log(\rho_{i})(\rho_{i}/\max_{j}\rho_{j})^{\alpha}}{\sum_{i=1}^{n}(\rho_{i}/\max_{j}\rho_{j})^{\alpha}}=-\log(\max_{i}\rho_{i}).$

This implies the Rényi entropy converges to the min-entropy.

(3) Tsallis entropy:

$\displaystyle
S_{q}(\rho)=\frac{1}{q-1}\left(1-\sum_{i=1}^{n}\rho_{i}^{q}\right),$

where $q\geq 1$. From the L’H$\hat{o}$spital theorem, we can also check that
Tsallis entropy converges to Shannon entropy as $q\searrow 1$.

From a simple calculation, we can also prove that both $H_{\alpha}(P)$ and
$S_{q}(P)$ is minimized when one of $p_{i}$ is one and the others are zero.
Since these properties are similar, we can expect that $\mathcal{F}$ can be
replaced by either Rényi entropy or Tsallis entropy. However, in this paper,
we focus on $\mathcal{F}$ defined in (3.2).

Now, we consider some entropies as $\mathcal{F}$ introduced in Example 4.1 as
potentials. Since we are considering two points graphs, we set
$\rho(r)=(r,1-r)$ with $0\leq r\leq 1$. Using the entropies that introduced in
Example 4.1, we define the functional $\mathcal{F}$ which can be expressed
(4.5) as follows:

(4.7)
$\displaystyle\mathcal{F}(r)=\begin{cases}-\Big{(}H(\rho(r))-H(\rho(1/2))\Big{)}=\log
2+r\log r+(1-r)\log(1-r)\quad&\text{Shannon entropy potential},\\\
-\Big{(}H_{\alpha}(\rho(r))-H_{\alpha}(\rho(1/2))\Big{)}=\log
2-\frac{1}{1-\alpha}\log(r^{\alpha}+(1-r)^{\alpha})\quad&\text{R\'{e}nyi
entropy potential},\\\
-\Big{(}S_{q}(\rho(r))-S_{q}(\rho(1/2))\Big{)}=\frac{1}{q-1}(r^{q}+(1-r)^{q}-2^{1-q})\quad&\text{Tsallis
entropy potential}.\end{cases}$

For the Shannon entropy potential, we subtract $H(\rho(1/2))$ from
$H(\rho(r))$. Since $\mathcal{F}(r)$ is positive and $H(\rho(r))-H(\rho(1/2))$
is negative, we put the minus sign. Not only the Shannon entropy, we also
applied this argument in Rényi entropy and Tsallis entropy to obtain Rényi
entropy potential and Tsallis entropy potential, respectively. From (4.5) and
(4.7), we can also define $\theta$, respectively. Relation (4.5) yields

$2\mathcal{F}(r)=\left(\int_{1/2}^{r}\frac{1}{\sqrt{\theta(r_{*})}}dr_{*}\right)^{2}.$

This implies that

$\int_{1/2}^{r}\frac{1}{\sqrt{\theta(r_{*})}}dr_{*}=\begin{cases}\displaystyle-\sqrt{2\mathcal{F}(r)}\quad&\displaystyle\text{if}\quad
0\leq r\leq\frac{1}{2},\\\
\displaystyle\sqrt{2\mathcal{F}(r)}\quad&\displaystyle\text{if}\quad\frac{1}{2}\leq
r\leq 1,\end{cases}$

and this again yields

$\frac{1}{\sqrt{\theta(r)}}=\begin{cases}\displaystyle-\frac{d}{dr}\sqrt{2\mathcal{F}(r)}\quad&\displaystyle\text{if}\quad
0\leq r\leq\frac{1}{2},\vspace{0.2cm}\\\
\displaystyle\frac{d}{dr}\sqrt{2\mathcal{F}(r)}\quad&\displaystyle\text{if}\quad\frac{1}{2}\leq
r\leq 1.\end{cases}$

Finally, we can obtain that

(4.8)
$\displaystyle\theta(r)=\left(\frac{d}{dr}\sqrt{2\mathcal{F}(r)}\right)^{-2}=\frac{2\mathcal{F}(r)}{\left(\frac{d}{dr}\mathcal{F}(r)\right)^{2}}.$

We substitute (4.7) into (4.8) to obtain

$\theta(r)=\begin{cases}\displaystyle\frac{2(\log 2+r\log
r+(1-r)\log(1-r))}{\left(\log r-\log(1-r)\right)^{2}}&\quad\text{Shannon
entropy},\\\
\displaystyle\frac{2(1-\alpha)^{2}(r^{\alpha}+(1-r)^{\alpha})^{2}\left(\log
2-\frac{1}{1-\alpha}\log(r^{\alpha}+(1-r)^{\alpha})\right)}{\alpha^{2}\left(r^{\alpha-1}-(1-r)^{\alpha-1}\right)^{2}}&\quad\text{R\'{e}nyi
entropy},\\\
\displaystyle\frac{2(q-1)(r^{q}+(1-r)^{q}-2^{1-q})}{q^{2}\left(r^{q-1}-(1-r)^{q-1}\right)^{2}}&\quad\text{Tsallis
entropy}.\end{cases}$

The minimum action $\mathcal{A}$ can be expressed as $\mathcal{F}$ as follows:

$\mathcal{A}(r_{0},r_{1})=\begin{cases}\displaystyle\frac{2\cosh 1}{\sinh
1}\left(\mathcal{F}(r_{0})+\mathcal{F}(r_{1})\right)+\frac{2}{\sinh
1}\sqrt{\mathcal{F}(r_{0})\mathcal{F}(r_{1})}\quad\text{if
}r_{0}\leq\frac{1}{2}\leq r_{1}~{}~{}\text{or}~{}~{}r_{1}\leq\frac{1}{2}\leq
r_{0},\vspace{0.2cm}\\\ \displaystyle\frac{2\cosh 1}{\sinh
1}\left(\mathcal{F}(r_{0})+\mathcal{F}(r_{1})\right)-\frac{2}{\sinh
1}\sqrt{\mathcal{F}(r_{0})\mathcal{F}(r_{1})}\quad\text{otherwise}.\end{cases}$

We notice that $\mathcal{A}$ is the value function for discrete Schrödinger
bridge system [6]. See also motivations in [3].

##  5\. Numerical examples

In this section, we provide some numeric examples. This section is consisted
of two parts. In the first part, we provide simulations of complete graphs. We
check (3.11) numerically. In the second part, we consider general graphs
instead of complete graphs. For this case, we have Proposition 3.1. We will
check this proposition again numerically on various discrete graphs.

###  5.1. Complete graph

In the first part of this subsection we fix parameters of (3.5) as follows:

$\Delta t=0.01,\quad\kappa=1,\quad n=4,\quad\rho^{0}=(0.5,0.3,0.15,0.05).$

We will check (3.11) numerically. First, we consider the case of
$\theta_{ij}=\min(\rho_{i},\rho_{j})$. Since $\rho_{1}^{0}$ is the unique
maximum of $\\{\rho_{1}^{0},\rho_{2}^{0},\rho_{3}^{0},\rho_{4}^{0}\\}$,
$\rho_{1}(t)$ converges to 1 exponentially. To check the exponential
convergence, we draw the graph of $\log(1-\rho_{1}(t))$. See Figure 1.

Figure 1. The case of $\theta_{ij}=\min(\rho_{i},\rho_{j})$. As we expected
$\log(1-\rho_{1}(t))$ converges to a line with a negative slope asymptotically
as $t\to\infty$.

Second, we consider the case of $\theta_{ij}=\min(\rho_{i},\rho_{j})^{2}$. We
know that $\rho_{1}(t)$ converges to 1 with the following convergence rate:
$1-\rho_{1}(t)\simeq\frac{1}{t}$. To check this convergence rate, we draw the
graph of $(1-\rho_{1}(t))^{-1}$. See Figure 2.

Figure 2. The case of $\theta_{ij}=\min(\rho_{i},\rho_{j})^{2}$. As we
expected the graph of $(1-\rho_{1}(t))^{-1}$ converges to a line with a
positive slope asymptotically as $t\to\infty$. This implies that
$\rho_{1}(t)\simeq\frac{1}{t}$.

Third, we consider the case of $\theta_{ij}=\min(\rho_{i},\rho_{j})^{3}$. We
know that $\rho_{1}(t)$ converges to 1 with the following convergence rate:
$1-\rho_{1}(t)\simeq\frac{1}{\sqrt{t}}$. To check this convergence rate, we
draw the graph of $(1-\rho_{1}(t))^{-2}$. See Figure 3.

Figure 3. The case of $\theta_{ij}=\min(\rho_{i},\rho_{j})^{3}$. As we
expected the graph of $(1-\rho_{1}(t))^{-2}$ converges to a line with a
positive slope asymptotically as $t\to\infty$. This implies that
$\rho_{1}(t)\simeq\frac{1}{\sqrt{t}}$.

###  5.2. General graph

In the second part of this subsection, we provide some numeric examples of
system (3.13) on general graphs. We consider three graphs containing six
edges($A$, $B$, $C$, $D$, $E$, and $F$). (1) Cycle graph(Figure 4), (2)
Lattice graph(Figure 5), and (3) Ribbon shaped graph(Figure 6). Each graph
contains six vertices and they have different topologies. In this part, we fix
the following parameters:

$\Delta
t=0.01,\quad\kappa=1,\quad\theta_{ij}=\min(\rho_{i},\rho_{j}),\quad\rho(0)=(0.3,0.2,0.1,0.1,0.1,0.2),$

where $\rho=(\rho_{A},\rho_{B},\rho_{C},\rho_{D},\rho_{E},\rho_{F})$.

ABCDEF Figure 4. Cycle graph.  ABCDEF Figure 5. Lattice graph CABDFE Figure 6.
Ribbon shaped graph

######  Example 5.1 (Cycle graph: Figure 4).

In this case, the density $\rho(t)$ converges to the following density:

$\rho^{\infty}=(0.7398,\quad 0,\quad 0,\quad 0.2602,\quad 0,\quad 0).$

We can observe that $A$ and $D$ are not connected, and the only nonzero part
of $\rho^{\infty}$ are $\rho_{A}^{\infty}$ and $\rho_{D}^{\infty}$. This
result follows Proposition 3.1.

######  Example 5.2 (Lattice graph: Figure 5).

If we remove one edge $AF$ from the cycle graph(Figure 4), then we can obtain
the lattice graph(Figure 5). In this case, the density $\rho(t)$ converges to
the following density:

$\rho^{\infty}=(0.5274,\quad 0,\quad 0,\quad 0.1958,\quad 0,\quad 0.2768).$

We can observe that $A$, $D$, and $F$ are not connected to each others, and
the only nonzero part of $\rho^{\infty}$ are $\rho_{A}^{\infty}$,
$\rho_{D}^{\infty}$, and $\rho_{F}^{\infty}$. This result follows Proposition
3.1.

######  Example 5.3 (Ribbon shaped graph: Figure 6).

In this case, the density $\rho(t)$ converges to the following density:

$\rho^{\infty}=(0.5948,\quad 0,\quad 0,\quad 0,\quad 0,\quad 0.4052).$

We can observe that $A$ and $F$ are not connected to each others, and the only
nonzero part of $\rho^{\infty}$ are $\rho_{A}^{\infty}$ and
$\rho_{F}^{\infty}$. This result follows Proposition 3.1.

###  5.3. Second order dynamics with $n\geq 3$.

In this subsection, we provide numeric results on second order dynamics with
$n\geq 3$. In this part, we fix the following parameters:

$\Delta t=0.01,\quad\kappa=1,\quad
n=6,\quad\theta_{ij}=\min(\rho_{i},\rho_{j})^{2}.$

We consider the initial data:

(5.1)
$\displaystyle\begin{cases}\rho^{0}=(0.3224,~{}0.2108,~{}0.1071,~{}0.0713,~{}0.2518,~{}0.0366),\\\
S^{0}=(0.1597,~{}-1.1129,~{}0.5929,~{}0.4568,~{}0.8299,~{}-0.2499).\end{cases}$

We present the numeric solution of system (3.3) with initial data (5.1) in
Figure 7.

Figure 7. Plot of $\rho_{i}(t)$(Left) and $S_{i}(t)$(Right) graphs with the
initial data (5.1). Since only one $\rho_{i}(t)$ converges to $1$ and others
converge to zero, we can observe the complete synchronization on the graph.

We also consider the other initial data:

(5.2)
$\displaystyle\begin{cases}\rho^{0}=(0.1524,~{}0.0910,~{}0.0698,~{}0.1583,~{}0.3424,~{}0.1862),\\\
S^{0}=(-0.4890,~{}-0.4542,~{}-0.2708,~{}-0.6929,~{}1.0627,~{}0.1228).\end{cases}$

We present the numeric solution of system (3.3) with initial data (5.2) in
Figure 8.

Figure 8. Plot of $\rho_{i}(t)$(Left) and $S_{i}(t)$(Right) graphs with the
initial data (5.2). Since only one $\rho_{i}(t)$ converges to $1$ and others
converge to zero, we can observe the complete synchronization on the graph.

## References

  * [1] Barb$\check{\mathrm{a}}$lat, I.: Systéms déquations différentielles d’oscillations non Linéaires, Rev. Math. Pures Appl. 4 (1959), 267-270.
  * [2] Chiarini, A., Conforti, G., Tamanini, L.: Schrödinger Problem for Lattice Gases: A Heuristic Point of View. In: Nielsen F., Barbaresco F. (eds) Geometric Science of Information. GSI 2021. Lecture Notes in Computer Science, Springer, Cham., 12829, (2021) 891-899.
  * [3] Chen, Y., Georgiou, T., and Pavon. M.: On the Relation Between Optimal Transport and Schrödinger Bridges: A Stochastic Control Viewpoint, J. Optim. Theory Appl., 169 (2016), 671-691.
  * [4] Chizat, L., Roussillon, P., Leger, F., Vialard, F.X., and Peyre, G.: Faster Wasserstein Distance Estimation with the Sinkhorn Divergence, NeurIPS (2020), 34th Conference on Neural Information Processing Systems.
  * [5] Chow, S.-N., Li, W., and Zhou, H.: A discrete Schrödinger equation via optimal transport on graphs, Journal of Functional Analysis. 276 (2019), 2440-2469.
  * [6] Chow, S.-N., Li, W., Mou, C., and Zhou, H.: Dynamical Schrödinger Bridge Problems on Graphs, Journal of Dynamics and Differential Equations (2021).
  * [7] Chow, S.-N., Li, W., and Zhou, H.: Entropy dissipation of Fokker-Planck equations on finite graphs, Discrete and Continuous Dynamical Systems-A, 38 (2018), 4929-4950.
  * [8] Chow, S.-N., Huang, W., Li, Y., and Zhou, H.: Fokker–Planck equations for a free energy functional or Markov process on a graph, Archive for Rational Mechanics and Analysis, 203 (2012), 969-1008.
  * [9] Ha, S.-Y, Ha, T., and Kim, J.-H.: On the complete synchronization of the Kuramoto phase model, Physics D, 239 (2010), 1692-1700.
  * [10] Ha, S.-Y., Kim, Y.-H., Morales, J., and Park, J. Y.: Emergence of phase concentration for the Kuramoto–Sakaguchi equation, Physica D: Nonlinear Phenomena. 401 (2020), 132154.
  * [11] Ha, S.-Y., Li, Z., and Xue, X.: Formation of phase-locked states in a population of locally interacting Kuramoto oscillators, Journal of Differential Equations, 255 (2013), 3053-3070.
  * [12] Kuramoto, Y.: Chemical Oscillations, Waves and Turbulence, Springer-Verlag, Berlin (1984)
  * [13] Kuramoto, Y.: Self-entrainment of a population of coupled non-linear oscillators, International symposium on mathematical problems in mathematical physics, Lecture notes in theoretical physics, 39 (1975), 420-422.
  * [14] Léger, F. and Li, W.: Hopf-Cole transformation via generalized Schrödinger bridge problem, Journal of Differential Equations. 274 (2021), 788-827.
  * [15] Li, W. and Montúfar, G.: Natural gradient via optimal transport, Information Geometry. 1 (2018), 181-214.
  * [16] Lohe, M. A.: Non-Abelian Kuramoto model and synchronization, J. Phys. A, 42 (2009), 395101.
  * [17] Mass, J.: Gradient flows of the entropy for finite Markov chains, Journal of Functional Analysis, 261 2250-2292, 2011.
  * [18] Mielke, A.: A gradient structure for reaction–diffusion systems and for energy-drift-diffusion, Nonlinearity, 24 (2011),13-29.
  * [19] Sakaguchi, H.: Cooperative phenomena in coupled oscillator systems under external fields, Prog. Theor. Phys. 79 (1988), 39-46.
  * [20] Villani, C.: Optimal transport: old and new, 338 (2008), Springer Science & Business media.
  * [21] Winfree, A.: Biological rhythms and the behavior of populations of coupled oscillators, J. Theoret. Bio., 16 (1967), 15-42.

##  Appendix A Detail Proofs

######  Lemma A.1.

Suppose that the initial data $(x^{0},p^{0})$ satisfy the following condition:

(A.1) $\displaystyle p^{0}=-\nabla U(x^{0}),$

and let $(x,p)$ be a solution to system (2.1). Then, we have

$p(t)=-\nabla U(x(t))\quad\forall~{}t\geq 0.$

###### Proof.

We have the following calculation:

$\displaystyle\frac{d}{dt}(p+\nabla U(x))_{i}$
$\displaystyle=\frac{d}{dt}p_{i}+\partial_{ij}U(x)\frac{dx_{j}}{dt}$
$\displaystyle=\frac{1}{2}\partial_{i}|\nabla
U(x)|^{2}+\partial_{ij}U(x)p_{j}$
$\displaystyle=\partial_{ij}U(x)\partial_{j}U(x)+\partial_{ij}U(x)p_{j}$
$\displaystyle=\partial_{ij}U(x)(p+\nabla U(x))_{j},$

where we use the Einstein’s convention. Then, we get the following inequality:

$\displaystyle\frac{1}{2}\frac{d}{dt}\|p+\nabla U(x)\|^{2}$
$\displaystyle=(p+\nabla U(x))_{i}\partial_{ij}U(x)(p+\nabla U(x))_{j}$
$\displaystyle\leq\|\nabla^{2}U(x)\|_{F}\cdot\|p+\nabla U(x)\|^{2}.$

From the initial condition $\|p(0)+\nabla U(x(0))\|=0$, we have

$\|p+\nabla U(x)\|^{2}=0$

for all $t\geq 0$. ∎

######  Lemma A.2.

Let $\rho_{t}$ be a solution to (2.5). If we set

(A.2) $\displaystyle
S_{t}=-\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}},$

then $(\rho_{t},S_{t})$ is a solution to (2.7) and the initial data satisfies
(2.8).

###### Proof.

Let $\rho_{t}$ be a solution to (2.5). From the substitution (A.2) and (2.5),
we get

$\partial_{t}\rho_{t}+\mathrm{div}(\rho_{t}\nabla S_{t})=0.$

So we have the first equation of (2.7). By direct calculations, we also have

$\displaystyle\partial_{t}S_{t}(x)$
$\displaystyle=-\int\frac{\delta^{2}\mathcal{F}(\rho_{t})}{\delta\rho_{t}^{2}}(x,y)\partial_{t}\rho_{t}(y)dy=-\int\frac{\delta^{2}\mathcal{F}(\rho_{t})}{\delta\rho_{t}^{2}}(x,y)\mathrm{div}\left(\rho_{t}\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)\right)(y)dy$
$\displaystyle=\int\rho_{t}(y)\frac{\delta}{\delta\rho_{t}(x)}\left|\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)(y)\right|^{2}dy$
$\displaystyle=\frac{\delta}{\delta\rho_{t}(x)}\left(\int\rho_{t}\left|\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)\right|^{2}dy\right)-\left|\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)(x)\right|^{2}$
$\displaystyle=\frac{\delta}{\delta\rho_{t}(x)}\left(\frac{1}{2}\int|\nabla
S_{t}|^{2}\rho_{t}dy\right)-\frac{1}{2}|\nabla S_{t}(x)|^{2}.$

Finally, we have

$\partial_{t}S_{t}+\frac{1}{2}|\nabla
S_{t}|^{2}=\frac{\delta}{\delta\rho_{t}}\left(\frac{1}{2}\int|\nabla
S_{t}|^{2}\rho_{t}dx\right).$

This implies that

$\partial_{t}S_{t}+\frac{1}{2}|\nabla
S_{t}|^{2}=\frac{\delta}{\delta\rho_{t}}\left(\frac{1}{2}\int\left|\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)\right|^{2}\rho_{t}(x)dx\right).$

Finally, the induced Wasserstein Hamiltonian flow from system (2.5) can be
written as follows:

$\displaystyle\begin{cases}\partial_{t}\rho_{t}+\mathrm{div}(\rho_{t}\nabla
S_{t})=0,\\\ \displaystyle\partial_{t}S_{t}+\frac{1}{2}|\nabla
S_{t}|^{2}=\frac{\delta}{\delta\rho_{t}}\left(\frac{1}{2}\int\left|\nabla\left(\frac{\delta\mathcal{F}(\rho_{t})}{\delta\rho_{t}}\right)\right|^{2}\rho_{t}(x)dx\right).\end{cases}$

This is the desired result. ∎

######  Lemma A.3.

Suppose that the initial condition $(\xi_{0},\xi^{*}_{0})$ satisfies

$\xi_{0}(x)\equiv 0\quad\forall~{}x\in M,$

and a pair of $C^{2}$ function $(\xi,\xi^{*})$ is a solution to system (2.10).
Then, $\xi_{t}(x)$ is identically zero for all $t\geq 0$ and $x\in M$.

###### Proof.

First, we define the following functional:

$\mathcal{A}(t)=\int\xi_{t}(x)^{2}dx.$

From the previous result (2.10), we have the following calculation:

$\displaystyle\partial_{t}\mathcal{A}$
$\displaystyle=2\int\xi_{t}(x)\nabla\xi_{t}(x)\cdot\nabla\xi^{*}(x)dx-\iint\xi_{t}(x)[\delta^{2}\mathcal{F}](x,u)\nabla_{u}\cdot(\rho_{t}(u)\nabla_{u}\xi_{t}(u))dxdu$
$\displaystyle=:\mathcal{I}_{1}+\mathcal{I}_{2}.$

Now, we estimate each terms for $0\leq t\leq T$. First, we estimate
$\mathcal{I}_{1}$:

$\displaystyle|\mathcal{I}_{1}|$
$\displaystyle=\left|\int\nabla\xi_{t}(x)^{2}\cdot\nabla\xi_{t}^{*}(x)dx\right|$
$\displaystyle=\left|\int\xi_{t}(x)^{2}\Delta\xi_{t}^{*}(x)dx\right|\leq\left(\sup_{t\in[0,T],x\in
M}|\Delta\xi_{t}^{*}(x)|\right)\mathcal{A}(t).$

Second, we estimate $\mathcal{I}_{2}$:

$\displaystyle|\mathcal{I}_{2}|$
$\displaystyle=\left|\iint\xi_{t}(x)\xi_{t}(u)\nabla_{u}\cdot(\rho_{t}(u)\nabla_{u}[\delta^{2}\mathcal{F}](x,u))dxdu\right|$
$\displaystyle\leq\left(\sup_{t\in[0,T],x,u\in
M}|\nabla_{u}\cdot(\rho_{t}(u)\nabla_{u}[\delta^{2}\mathcal{F}](x,u))|\right)\mathcal{A}(t)\mathrm{Vol}(M).$

From the above estimates, for any $T>0$, we can find a fixed positive constant
$C_{T}$ which satisfies the follows:

$\partial_{t}\mathcal{A}(t)\leq C_{T}\mathcal{A}(t)\quad\forall 0\leq t\leq
T.$

This property yields

$\mathcal{A}(t)\leq\mathcal{A}(0)e^{C_{T}t}\quad\forall~{}0\leq t\leq T.$

From the initial condition $\mathcal{A}(0)=0$, we get $\mathcal{A}(t)=0$ for
all $0\leq t\leq T$. Since this argument holds for all $T>0$, we can prove
that $\mathcal{A}(t)=0$ for all $t\geq 0$. This yields $\xi_{t}(x)\equiv 0$
for all $x\in M$ and $t\geq 0$. ∎