# Phase transition in a kinetic mean-field game model of inertial self-
propelled agents

Piyush Grover<EMAIL_ADDRESS>University of Nebraska-Lincoln  Mandy Huo
California Institute of Technology

###### Abstract

The framework of Mean-field Games (MFGs) is used for modelling the collective
dynamics of large populations of non-cooperative decision-making agents. We
formulate and analyze a kinetic MFG model for an interacting system of non-
cooperative motile agents with inertial dynamics and finite-range
interactions, where each agent is minimizing a biologically inspired cost
function. By analyzing the associated coupled forward-backward in time system
of nonlinear Fokker-Planck and Hamilton-Jacobi-Bellman equations, we obtain
conditions for closed-loop linear stability of the spatially homogeneous MFG
equilibrium that corresponds to an ordered state with non-zero mean speed.
Using a combination of analysis and numerical simulations, we show that when
energetic cost of control is reduced below a critical value, this equilibrium
loses stability, and the system transitions to a travelling wave solution. Our
work provides a game-theoretic perspective to the problem of collective motion
in non-equilibrium biological and bio-inspired systems.

> How can we understand phase transitions in collective motion of a population
> of motile agents making individually optimal decisions over time ? Mean
> field game (MFG) theory is a mathematical framework for modelling such
> systems. We study a system consisting of a large number of agents that
> optimally modulate their acceleration to minimize the long term average of a
> weighted sum of two costs. The first is the total energy use, and the second
> penalizes the mismatch between the local mean speed and a preferred speed.
> We show that the system transitions from a spatially homogeneous ordered
> state to a travelling wave as the energetic cost is reduced below a critical
> value.

##  I Introduction

This paper is concerned with phase transitions in the collective motion of
optimal decision-making non-cooperative motile agents. Collective motion and
synchronization continues to be actively studied in many scientific
domainsVicsek  _et al._ (1995); Levine, Rappel, and Cohen (2000); Strogatz
(2000); Cucker and Smale (2007); Gompper  _et al._ (2020); Kassabov, Strogatz,
and Townsend (2021). In purely physical systems, the interactions between
particles are given by physical laws, and systems are studied using agent-
based (or microscopic), kinetic and hydrodynamics models that incorporate
those laws. On the other hand, models for systems involving decision-making
_agents_ , e.g., bird flocking, traffic, human crowds and robot swarms, often
incorporate other types of interactions such as collision avoidance,
alignment, and cohesion towards the mean. Most commonly, these latter types of
interactions are captured in phenomenological models using physical analogies.

A different, ‘inverse modelling’ approach Ding  _et al._ (2022) is to
stipulate that the collective behavior of a population of decision-making
agents is a solution to a collective optimization or optimal control problem.
When considering large number of non-cooperative agents making sequential
decisions, the framework of Mean-Field Games (MFGs) Caines (2021) is
appropriate. In a MFG system, the collective behavior is the result of each
agent solving an optimal control problem that depends on its own state and
control as well as the collective state Bensoussan, Frehse, and Yam (2013);
Huang, Caines, and Malhame (2007). MFGs formulated in continuous state space
and time are described by coupled set of forward-backward in time nonlinear
partial differential equations (PDEs). While standard kinetic or hydrodynamic
equations used for modelling collective behavior are initial value problems
(IVP or evolution PDEs) , the MFG systems have a forward-backward in time
structure, and hence consist of boundary value problem (BVP in time PDEs) .

In evolution PDEs, phase transitions Solon  _et al._ (2015); Strogatz (2000);
Orosz, Krauskopf, and Wilson (2005) in collective behavior are studied via
stability and bifurcation analysis of steady or time-periodic states. In
recent years, similar questions are being explored in the context of MFGs,
where bifurcations can be interpreted as phase transitions in the collective
behavior of optimal decision-making agents as problem parameters such as the
unit cost of control are varied. While uniqueness of solutions in MFGs is
guaranteed under certain monotonicity conditions Bardi and Fischer (2019),
there is an increasing interest Nourian, Caines, and Malhamé (2011); Grover,
Bakshi, and Theodorou (2018); Cirant (2019); Cirant, Nurbekyan  _et al._
(2018); Ullmo, Swiecicki, and Gobron (2019); Salhab, Malhame, and Le Ny
(2019); Hongler (2020); Lori and Grover (2024) in studying MFGs that don’t
satisfy such conditions and hence are expected to have multiple co-existing
steady states.

In this context, starting with the seminal work of Yin et al. Yin  _et al._
(2012), several authors have studied the phase transition from a ‘disordered’
steady state to an ‘ordered’ travelling wave solution in the MFG problem
involving synchronization of oscillatorsCarmona, Cormier, and Soner (2023);
Cesaroni and Cirant (2024). The system involves globally cost-coupled agents
with first-order dynamics (i.e., without inertia), moving on a circular
domain. Hence, these MFG models can be understood as ‘inverse models’ for the
Kuramoto dynamics Kuramoto (1975). Other works Cirant (2019); Cirant,
Nurbekyan  _et al._ (2018) have proved the existence of time-periodic
solutions in certain classes of MFGs set in periodic spatial domains.

In this paper, we study phase transitions in a MFG ‘inverse model’ for the
kinetic Czirók model Czirók, Barabási, and Vicsek (1999), where unlike all the
works discussed above, the population consists of locally interacting inertial
agents, i.e., the agent dynamics are _second-order in time_ , and the agents
control their acceleration. As a result, the ordered states of the Czirók
model include spatially homogeneous equilibria and spatially inhomogeneous
travelling waves. We show that the MFG model we derive inherits these
features. Specifically, we provide a stability analysis of the spatially
homogeneous equilibria of the kinetic MFG Degond, Liu, and Ringhofer (2014).
Further, we numerically show the phase transition in the MFG to travelling
waves as the equilibrium states lose stability when the unit cost of control
is reduced below a critical value. Our analysis combines techniques use to
study generalized Kuramoto models with inertia Barre and Métivier (2016);
Acebrón, Bonilla, and Spigler (2000), and invariant subspace methods employed
in the study of Riccati equations Laub (1991); Abou-Kandil  _et al._ (2012) as
well as linear-quadratic (LQ) MFGs Huang and Zhou (2019).

##  II Czirók model revisited

In this section, we briefly review the Czirók Garnier, Papanicolaou, and Yang
(2019) model, and provide a new stability analysis which sets the stage for
the stability analysis of the MFG model in the later sections. The model is a
biologically inspired interacting particle system of $n$ agents, where each
agent is moving on a 1D periodic domain $[0,l]$. The second-order dynamics of
the $i$th agent are given by the stochastic differential equation (SDE):

$\displaystyle\textrm{d}x_{i}$ $\displaystyle=u_{i}\textrm{d}t,$ (1)
$\displaystyle\textrm{d}u_{i}$ $\displaystyle=[G(\langle
u\rangle_{i})-u_{i}]\textrm{d}t+\sigma\textrm{d}w_{i}(t),$ (2)

where $\langle
u\rangle_{i}=\dfrac{1}{n}\displaystyle\sum^{n}_{j=1}u_{j}\phi(\|x_{j}-x_{i}\|)$,
and $\|x_{j}-x_{i}\|=\min(|x_{j}-x_{i}|,l-|x_{j}-x_{i}|)$ is the distance on
the torus between the positions of the $i^{\textrm{th}}$ and $j$th agents. The
$w_{i}$ are independent Brownian motions, and $\sigma$ is noise intensity. The
finite-range interaction kernel $\phi(x)=\dfrac{l}{2}\times\mathbf{1}_{[0,1]}$
with $\frac{1}{l}\bigintssss^{l}_{0}\phi(\|x\|)\textrm{d}x=1$, and
$G(u)=\dfrac{h+1}{5}u-\dfrac{h}{125}u^{3}$. Here, $G(u)-u$ is the gradient of
a bistable potential whose critical points are the preferred mean speeds. In
the $n\rightarrow\infty$ limit, the nonlinear PDE describing the evolution of
density $\rho(t,x,u)$ is given by the Fokker-Planck equation

$\displaystyle\frac{\partial{\rho(t,x,u)}}{\partial t}=$
$\displaystyle\\-u\frac{\partial{\rho(t,x,u)}}{\partial
x}-\frac{\partial}{\partial
u}\left(\left[G\left(\int_{0}^{L}\int_{-\infty}^{\infty}u^{\prime}\phi(\|x^{\prime}\|)\rho(t,x-x^{\prime},u^{\prime})du^{\prime}dx^{\prime}\right)-u\right]\rho(t,x,u)\right)$
(3)
$\displaystyle+\frac{1}{2}\sigma^{2}\frac{\partial^{2}\rho(t,x,u)}{\partial
u^{2}}.$

Note that this equation is degenerate since there is no diffusion in the
spatial ($x$) variable. The system has spatially homogeneous equilibrium
states $\rho_{\xi}(u)$ of the form $\rho_{\xi}(u)=\dfrac{1}{l}F_{\xi}(u)$,
where
$F_{\xi}(u)=\dfrac{1}{\sqrt{\pi\sigma^{2}}}e^{-(u-\xi)^{2}/{\sigma^{2}}}$, and
$\xi$ is the mean speed of the population satisfying the fixed point equation
$\xi=G(\xi)$. This fixed point equation has a unique ‘disordered’ equilibrium
with $\xi=0$ for $h<4$, and two additional ordered equilibria with $\pm\xi\neq
0$ for $h>4$.

###  II.1 Linearization and the operator eigenvalue equation

We focus our stability analysis on the ordered ($\xi\neq 0$) spatially
homogeneous equilibria which exist when $h>4$. It is known that for $h>4$, the
disordered equilibrium is unstable. We linearize Eq. 3 around $\rho_{\xi}(u)$,
and note that $F^{\prime}_{\xi}(u)=-\dfrac{2(u-\xi)}{\sigma^{2}}F_{\xi}(u)$.
We take the perturbed state to be of the form

$\displaystyle\rho(t,x,u)=\rho_{\xi}(u)+\epsilon\sqrt{\rho_{\xi}(u)}\bar{\rho}(t,x,u),$
(4)

and expand the perturbation $\bar{\rho}$ using spatial Fourier decomposition:

$\displaystyle\bar{\rho}(t,x,u)=\sum_{k=0,\pm 1,\pm
2,...}\hat{\rho}_{k}(t,u)e^{i2\pi kx/l}.$ (5)

We will work in the space of square-integrable complex-valued functions on the
real line, i.e, we will assume $\hat{\rho}(t,.)\in L^{2}(\mathbb{R},du)$, with
inner product $\langle
a,b\rangle=\bigintssss_{\infty}^{\infty}\bar{a}(u)b(u)du$. The resulting
linearized evolution equation for the $k$th Fourier mode is

$\displaystyle\frac{\partial{\hat{\rho}_{k}(t,u)}}{\partial
t}=L^{k}\hat{\rho}_{k}=(L^{k}_{loc,1}+L^{k}_{nonloc})\hat{\rho}_{k},$ (6)

where

$\displaystyle L^{k}_{loc,1}[f](u)=\left[\frac{-i2\pi
ku}{l}-\frac{[(u-\xi)^{2}-\sigma^{2}]}{2\sigma^{2}}\right]f(u)+\frac{1}{2}\sigma^{2}\frac{\partial^{2}f(u)}{\partial
u^{2}},$ (7) $\displaystyle
L^{k}_{nonloc}[f](u)=\frac{2(u-\xi)}{\sigma^{2}}\sqrt{F_{\xi}(u)}G^{\prime}(\xi)\phi_{k}\left(\int
u^{\prime}\sqrt{F_{\xi}(u^{\prime})}f(u^{\prime})du^{\prime}\right).$ (8)

The eigenfunctions and eigenvalues of the operator $L^{k}_{loc,1}$ can be
computed explicitly as follows. We perform a change of variables Acebrón,
Bonilla, and Spigler (2000) $u=\mathfrak{g}(v)=a_{1}v+a_{2}$, where
$a_{1}=\sigma/\sqrt{2}$ , and $a_{2}=\xi-2ik\pi\sigma^{2}/l$. With this
substitution, the eigenvalue equation for $L^{k}_{loc,1}$ reduces to

$\displaystyle\alpha_{k}\hat{\rho}_{loc,k}=\left[\left(\frac{1}{2}-\frac{v^{2}}{4}+c_{2}(k)\right)\hat{\rho}_{loc,k}+\frac{\partial^{2}\hat{\rho}_{loc,k}}{\partial
v^{2}}\right].$ (9)

It is known Abramowitz and Stegun (1964) that the eigenvalue Eq. 9 has
solutions $\alpha_{k,p}=-p+c_{2}(k),p=0,1,2,\dots$. Here
$c_{2}(k)=-\dfrac{2\pi^{2}k^{2}\sigma^{2}}{l^{2}}-\dfrac{2ik\pi\xi}{l}.$ The
corresponding eigenfunctions are parabolic cylinder functions
$D_{p}(v)=2^{-\frac{p}{2}}e^{\frac{-v^{2}}{4}}H_{p}(\dfrac{v}{\sqrt{2}})=e^{\frac{-v^{2}}{4}}\tilde{H}_{p}(v)$.
Here $H_{p}$ are the physicist’s Hermite polynomials and $\tilde{H}_{p}$ are
the probabilist’s Hermite polynomials. Hence, operator $L^{k}_{loc,1}$ also
has the same eigenvalues $\alpha_{k,p}$. Its eigenfunctions are
$\eta_{k,p}(u)=z(p)D_{p}(\mathfrak{g}^{-1}(u))=z(p)D_{p}(\dfrac{\sqrt{2}}{\sigma}(u-\xi+\dfrac{2ik\pi\sigma^{2}}{l}))$.
Here $z(p)=\sqrt{\dfrac{1}{\sqrt{2\pi}a_{1}p!}}$ is a normalizing factor. The
adjoint $(L^{k}_{loc,1})^{*}$ has eigenvalues $\bar{\alpha}_{k,p}$, and
eigenfunctions $\psi_{k,p}=\bar{\eta}_{k,p}$. Note that due to our chosen
normalization, $\langle\eta_{k,p},\psi_{k,q}\rangle=\delta_{pq}$. Hence the
set $\\{\eta_{k,p},\psi_{k,q}\\}_{p=0,1,\dots,q=0,1,\dots}$ form a
biorthogonal eigenfunction basis set of $L^{2}[\mathbb{R},du]$.

With the eigenfunctions and eigenvalues of the local operator at hand, we can
derive the characteristic equation of full operator $L^{k}$ as follows. Let
($\lambda,\hat{\rho}$) be an eigenvalue-eigenfunction pair of $L^{k}$. The
eigenvalue equation is

$\displaystyle\lambda\hat{\rho}(u)=L^{k}_{loc,1}[\hat{\rho}](u)+g_{0}(u)\langle\bar{\hat{\rho}},s\rangle,$
(10)

where we have defined functions
$g_{0}(u)\triangleq\dfrac{2(u-\xi)}{\sigma^{2}}\sqrt{F_{\xi}(u)}G^{\prime}(\xi)\phi_{k}$,
and $s(u)\triangleq u\sqrt{F_{\xi}(u)}$.

Let $R^{k}_{1,\lambda}=(L^{k}_{loc,1}-\lambda)^{-1}$ be the resolvent of
$L^{k}_{loc,1}$. Then, the action of $R^{k}_{1,\lambda}$ on an arbitrary
complex-valued function f(u) is
$R^{k}_{1,\lambda}[f](u)=\displaystyle\sum_{p=0}^{\infty}\dfrac{\langle\psi_{k,p},f\rangle}{\alpha_{k,p}-\lambda}\eta_{k,p}(u)$.
Using this relation, Eq. 10 can be rewritten as
$\hat{\rho}=-\displaystyle\sum_{p=0}^{\infty}\dfrac{\langle\psi_{k,p},g_{0}\rangle\langle\bar{\hat{\rho}},s\rangle}{\alpha_{k,p}-\lambda}\eta_{k,p}(u)$.
Taking the inner product of the conjugate of this expression with $s$,
cancelling common terms on both sides, and another conjugate operation yields
the characteristic equation for $L^{k}_{loc,1}$:

$\displaystyle
1=\sum_{p=0}^{\infty}\frac{\langle\psi_{k,p},g_{0}\rangle\langle\bar{\eta}_{k,p},s\rangle}{\lambda-\alpha_{k,p}}=\sum_{p=0}^{\infty}\frac{\langle\bar{\eta}_{k,p},g_{0}\rangle\langle\bar{\eta}_{k,p},s\rangle}{\lambda-\alpha_{k,p}}.$
(11)

To verify our analysis, we solve Eq. 11 for $l=10,h=5$, and $0.2\leq\sigma\leq
2.5$ using 20 eigenfunctions for each Fourier mode. As shown in Fig. 1, the
$k=1$ Fourier mode loses stability at $\sigma=\sigma_{c}=1.8$, which matches
the prediction in Garnier et al. Garnier, Papanicolaou, and Yang (2019).

(a)

(b)

Figure 1:  (a) The real (solid) and imaginary (dashed) components of the
eigenvalue of $L_{k}$ (for $h=5,k=1$) closest to the imaginary axis, as a
function of noise intensity $\sigma$. The system is unstable for
$\sigma\leq\sigma_{c}=1.8$ as the real part is positive in that range. (b) The
spectrum of $L_{k}$ at the threshold of stability, $\sigma=\sigma_{c}$.

##### Transition to travelling wave:

When a spatially homogeneous equilibrium $\rho_{\xi}(u)$ ($\xi\neq 0$) loses
stability, the system converges to a stable travelling wave Garnier,
Papanicolaou, and Yang (2019), i.e., a time-varying spatially non-homogeneous
solution of the form $\rho(t,x,v)=\tilde{\rho}(x-\omega t,v)$. In Fig. 2, we
show one such solution obtained by numerically solving the IVP Eqs. 3 in the
open-source software Dedalus Burns  _et al._ (2020).

Figure 2: The marginal density $\bigintssss_{\mathbb{R}}\rho(t,x,v)dv$ of a
travelling wave solution of the Czirók model Eq. 3 for
$h=5,\sigma=0.8<\sigma_{c}$. This solution is the steady state reached upon
perturbing the unstable spatially homogeneous equilibrium $\rho_{\xi}(u)$,
where $\bigintssss_{\mathbb{R}}\rho_{\xi}(v)dv=\dfrac{1}{l}=0.1.$

##  III A Mean-field game Czirók model

To construct the kinetic MFG model, we consider the interacting particle
system of $n$ agents, in which the $i$th agent has the dynamics of the form

$\displaystyle\textrm{d}x_{i}$ $\displaystyle=u_{i}\textrm{d}t,$ (12)
$\displaystyle\textrm{d}u_{i}$
$\displaystyle=\alpha_{i}\textrm{d}t+\sigma\textrm{d}w_{i}(t),$ (13)

where $(x_{i}(t),u_{i}(t))$ are its position and speed, and $\alpha_{i}(t)$ is
its chosen control. Each agent seeks to minimize the long term average of the
sum of two costs. The first cost depends on a weighted mean speed of the
population around the agent’s position, and the other is a measure of the
control effort. The cost function inspired by the interaction term in the
Czirók model Eq. (3) is

$J=\limsup_{T\rightarrow\infty}\frac{1}{T}\int^{T}_{0}[c((x_{i},u_{i});(x,u)_{-i})+r\alpha_{i}^{2}]\textrm{d}t,$
(14)

with cost-coupling function

$c((x_{i},u_{i});(x,u)_{-i})=\left([G(\langle
u\rangle_{i})-u_{i}]\right)^{2},$ (15)

where
$(x,u)_{-i}=\\{(x_{1},u_{1}),(x_{2},u_{2}),\dots,(x_{i-1},u_{i-1}),(x_{i+1},u_{i+1}),\dots,(x_{n},u_{n})\\}$
represents the rest of the population. Here, $r$ is the unit cost of control.
In the $n\rightarrow\infty$ population limit, we can approximate the cost-
coupling function as Grover, Bakshi, and Theodorou (2018)

$c[\rho](t,x,u)=\left[G\left(\iint
u^{\prime}\phi(\|x^{\prime}\|)\rho(t,x-x^{\prime},u^{\prime})\textrm{d}u^{\prime}\textrm{d}x^{\prime}\right)-u\right]^{2}.$
(16)

Using standard techniques Yin  _et al._ (2012); Nourian (2013), we can show
that the corresponding MFG system consist of the following nonlinear Fokker-
Planck (FP) and Hamilton-Jacobi-Bellman (HJB) equations governing the density
$\rho(x,u,t)$ and the relative value function $h(x,u,t)$ respectively:

$\displaystyle\frac{\partial{\rho}}{\partial
t}=-u\frac{\partial{\rho}}{\partial x}+\frac{1}{2r}\frac{\partial}{\partial
u}\\{\rho\frac{\partial h}{\partial
u}\\}+\frac{\sigma^{2}}{2}\frac{\partial^{2}\rho}{\partial u^{2}},$ (17)
$\displaystyle\frac{\partial{h}}{\partial t}=\chi-c[\rho]-u\frac{\partial
h}{\partial x}+\frac{1}{4r}(\frac{\partial h}{\partial
u})^{2}-\frac{\sigma^{2}}{2}\frac{\partial^{2}h}{\partial u^{2}},$ (18)

where $\chi$ is the minimum average cost. The PDEs can be solved given an
initial density $\rho(0.,.)$, the required final condition
$\displaystyle\lim_{t\rightarrow\infty}h(t,.,.,)=0$, and appropriate decay
conditions in space. The optimal control
$\alpha(t,x,u)=-\dfrac{1}{2r}\dfrac{\partial h}{\partial u}$. It is easily
verified that the $(\rho_{\xi^{*}}(u),h_{\xi^{*}}(u))$ is a stationary
solution to the MFG equations, where
$\rho_{\xi^{*}}(u)=\dfrac{1}{l}F_{\xi^{*}}(u),F_{\xi^{*}}(u)=\dfrac{1}{\sqrt{\pi\sqrt{r}\sigma^{2}}}e^{-(u-\xi^{*})^{2}/{(\sqrt{r}\sigma^{2})}}$
and $h_{\xi^{*}}(u)=\sqrt{r}(u-\xi^{*})^{2}$. Here $\xi^{*}$ is again the
mean-speed of the population, satisfying $G(\xi^{*})=\xi^{*}$, and
$\chi=\sigma^{2}\sqrt{r}$. Hence, similar to Czirók model discussed in Sec.
II, the MFG model also possesses three spatially homogeneous equilibria for
$h>4$, among which the two equilibria with $\pm\xi^{*}\neq 0$ are the ordered
states.

###  III.1 Linearization and the operator eigenvalue equation

We linearize the MFG Eqs. (17,18) around an ordered equilibrium
$(\rho_{\xi}(u),h_{\xi}(u))$ with $\xi=\xi^{*}\neq 0$. As in section II.1, we
restrict the stability analysis to a class of density perturbations that are
exponentially decaying in $u$, i.e., we choose
$\rho(t,x,u)=\rho_{\xi}(u)+\epsilon\sqrt{\rho_{\xi}(u)}\bar{\rho}(t,x,u)$,
where $\bar{\rho}(t,x,.)\in\mathbb{L}^{2}\\{\mathbb{R},du\\}$. In case of the
value function, the perturbations are restricted to have a bounded exponential
growth rate, i.e., we choose
$h(t,x,u)=h_{\xi}(u)+\epsilon\dfrac{\bar{h}(t,x,u)}{\sqrt{\rho_{\xi}(u)}}$,
where $\bar{h}(t,x,.)\in\mathbb{L}^{2}\\{\mathbb{R},du\\}.$ We expand the
perturbations using spatial Fourier decompositions,
$\bar{\rho}(t,x,u)=\displaystyle\sum_{k=0,\pm 1,\pm
2,...}\hat{\rho}_{k}(t,u)e^{i2\pi kx/l}$, and
$\bar{h}(t,x,u)=\displaystyle\sum_{k=0,\pm 1,\pm
2,...}\hat{h}_{k}(t,u)e^{i2\pi kx/l}$. With these substitutions, the
linearized FP equation is

$\displaystyle\frac{\partial{\hat{\rho}_{k}(t,u)}}{\partial t}=-\frac{i2\pi
ku}{l}\hat{\rho}_{k}-\tilde{L}_{loc}\hat{\rho}_{k}-\frac{1}{r\sigma^{2}}\tilde{L}_{loc}\hat{h}_{k},$
(19)

where
$\tilde{L}_{loc}\triangleq\left[\dfrac{[(u-\xi)^{2}-\sqrt{r}\sigma^{2}]}{2r\sigma^{2}}\right]-\dfrac{\sigma^{2}}{2}\dfrac{\partial^{2}}{\partial
u^{2}}$ is a self-adjoint operator. Hence, $L^{k}_{loc,1}=-\dfrac{i2\pi
ku}{l}-\tilde{L}_{loc}$.

The linearized HJB equation is

$\displaystyle\frac{\partial{\hat{h}}_{k}(t,u)}{\partial t}=-\frac{i2\pi
ku}{l}\hat{h}_{k}(t,u)+\tilde{L}_{loc}\hat{h}_{k}(t,u)-\tilde{c}[\hat{\rho}_{k}],$
(20)

where
$\tilde{c}[\hat{\rho}_{k}]=2(\xi-u)G^{\prime}(\xi)\phi_{k}\sqrt{F_{\xi}(u)}\left(\bigintssss
u^{\prime}\sqrt{F_{\xi}(u^{\prime})}\hat{\rho}_{k}(t,u^{\prime})du^{\prime}\right)$.
We define $L^{k}_{loc,2}\triangleq-\dfrac{i2\pi
ku}{l}+\tilde{L}_{loc}=-(L^{k}_{loc,1})^{*}$. The eigenvalues of
$L^{k}_{loc,2}$ are
$\beta_{k,p}=-\bar{\alpha}_{k,p}=\dfrac{p}{\sqrt{r}}-\dfrac{i2\pi
k\xi}{l}+\dfrac{2k^{2}\pi^{2}\sigma^{2}r}{l^{2}}$. The corresponding
eigenfunctions are $\psi_{k,p}(u)$. The coupled system of equations for the
linearized MFG are

$\displaystyle\frac{\partial{\hat{\rho}}_{k}(t,u)}{\partial
t}=L^{k}_{loc,1}{\hat{\rho}}_{k}(t,u)-\frac{1}{r\sigma^{2}}\tilde{L}_{loc}{\hat{h}}_{k}(t,u),$
(21) $\displaystyle\frac{\partial{\hat{h}}_{k}(t,u)}{\partial
t}=-\tilde{c}[\hat{\rho}_{k}]+L^{k}_{loc,2}{\hat{h}}_{k}(t,u).$ (22)

Let $(\\{\hat{\rho}_{k},\hat{h}_{k}\\},\lambda)$ be a eigenfunction-eigenvalue
pair for the linearized MFG system. The abstract operator eigenvalue equation
can be written as:

$\displaystyle\lambda\begin{bmatrix}\hat{\rho}_{k}\\\ \\\
\hat{h}_{k}\end{bmatrix}=\begin{bmatrix}L^{k}_{loc,1}&-\frac{1}{r\sigma^{2}}\tilde{L}_{loc}\\\
\\\ 0&L^{k}_{loc,2}\end{bmatrix}\begin{bmatrix}\hat{\rho}_{k}\\\ \\\
\hat{h}_{k}\end{bmatrix}+\begin{bmatrix}0\\\ \\\
-g(u)\langle\bar{\hat{\rho}}_{k},s\rangle\end{bmatrix},$ (23)

where, $g(u)\triangleq 2(\xi-u)\sqrt{F_{\xi}(u)}G^{\prime}(\xi)\phi_{k}$, and
$s(u)=u\sqrt{F_{\xi}(u)}$. Analogous to the development in Sec. II.1, we
employ the resolvents $R^{k}_{1,\lambda}$ of $L^{k}_{loc,1}$, and
$R^{k}_{2,\lambda}$ of $L^{k}_{loc,2}$ to obtain the characteristic equation
for the MFG eigenvalues. Here,
$R^{k}_{2,\lambda}[f](u)=\displaystyle\sum_{q=0}^{\infty}\dfrac{\langle\eta_{q},f\rangle}{\beta_{k,q}-\lambda}\psi_{k,q}(u)$.
Using this relation, the second equation of Eqs. 23 can be written as
$\hat{h}_{k}=\displaystyle\sum_{q=0}^{\infty}\dfrac{\langle\eta_{q},g\rangle\langle\bar{\hat{\rho}}_{k},s\rangle}{\beta_{k,q}-\lambda}\psi_{k,q}(u)$.
Substituting this expression of $\hat{h}_{k}$ in the first equation of Eqs. 23
yields

$\displaystyle
R^{k}_{1,\lambda}[\hat{\rho}_{k}](v)=\dfrac{1}{r\sigma^{2}}\sum_{q=0}^{\infty}\dfrac{\langle\eta_{q},g\rangle\langle\bar{\hat{\rho}}_{k},s\rangle}{\beta_{k,q}-\lambda}\tilde{L}_{loc}[\psi_{k,q}](u).$
(24)

After applying $(R^{k}_{1,\lambda})^{-1}$ to above equation, we take the inner
product of the conjugate of the resulting equation with $s$. Cancelling common
factors from both sides, and taking the conjugate again yields the
characteristic equation for the MFG eigenvalues:

$\displaystyle
1=\frac{1}{r\sigma^{2}}\sum_{p=0}^{\infty}\sum_{q=0}^{\infty}\frac{\langle\bar{\eta}_{k,p},s\rangle\langle\bar{\psi}_{q},g\rangle\langle\bar{\eta}_{k,p},\tilde{L}_{loc}[\psi_{k,q}]\rangle}{(\alpha_{k,p}-\lambda)(\beta_{k,q}-\lambda)}.$
(25)

###  III.2 Linear stability of the ordered MFG equilibria

Since a MFG system is a forward-backward system, the stability of an
equilibrium (or more general invariant sets such as periodic orbits or
travelling waves) needs to be defined in a way that captures the feedback
nature of the system. Consider the linear two-point BVP on the time interval
$[0,\infty]$ consisting of the linearized MFG Eqs. (21,22) along with a
prescribed initial density perturbation $\hat{\rho}_{k}(0,.)$, and final value
function $\displaystyle\lim_{t\rightarrow\infty}\hat{h}_{k}(t,.)=0$ (inherited
from the HJB equation).

######  Definition 1.

Huang, Malhamé, and Caines (2006); Guéant (2009); Yin  _et al._ (2012);
Grover, Bakshi, and Theodorou (2018) A MFG equilibrium $(\rho_{\xi},h_{\xi})$
of the MFG system (17,18) is said to be linearly stable if the above defined
BVP has a unique solution $(\hat{\rho}_{k}(t,.),\hat{h}_{k}(t,.))$ such that
$\displaystyle\lim_{t\rightarrow\infty}\hat{\rho}_{k}(t,.)\rightarrow 0$. In
other words, any initial density perturbation decays to $0$ in the closed loop
as the time horizon goes to $\infty$.

To proceed with the stability analysis, we first obtain the explicit form of
the linear BVP discussed above, using the eigenfunction expansions
$\hat{\rho}_{k}(t,u)=\displaystyle\sum_{p}Y_{1,p}(t)\eta_{k,p}(u)$, and
$\hat{h}_{k}(t,u)=\displaystyle\sum_{p}Y_{2,p}(t)\psi_{k,p}(u)$. The
prescribed initial condition on density fixes the value of $Y_{1}(0)$, and the
final condition on the value function becomes
$\displaystyle\lim_{t\rightarrow\infty}Y_{2}(t)=0$. The ODE system obtained by
inserting the expansions into Eqs. (21,22) is

$\displaystyle\begin{bmatrix}\dot{Y}_{1}\\\
\dot{Y}_{2}\end{bmatrix}=N\begin{bmatrix}Y_{1}\\\ Y_{2}\end{bmatrix},$ (26)

where $N\triangleq\begin{bmatrix}A_{1}&B_{1}\\\ A_{2}&B_{2}\end{bmatrix}$. The
submatrices are $(A_{1})_{q,p}=\alpha_{k,q-1}\delta_{qp}$,
$(A_{2})_{q,p}=-\langle\eta_{k,q-1},g(u)\rangle\langle
s(u),\eta_{k,p-1}\rangle$,
$(B_{1})_{q,p}=\dfrac{-1}{r\sigma^{2}}\langle\psi_{k,q-1},\tilde{L}_{loc}\psi_{k,p-1}\rangle$,
and $(B_{2})_{q,p}=\beta_{k,q-1}\delta_{qp}$. Using the properties of
operators involved, we make the following observations. First, $A_{1}$ is
diagonal, $Re(\sigma(A_{1}))<0$, i.e., all eigenvalues of $A_{1}$ lie in the
left half plane, and $B_{2}=-A_{1}^{*}$. Second, $B_{1}$ is Hermitian, hence
its eigenvalues are real, with $Re(\sigma(B_{1}))\leq 0$. Finally, $A_{2}$ is
a rank-1 non-Hermitian matrix. The characteristic equation of
$N=\Delta+\mathcal{P}_{1}$ is

$\displaystyle
det(\Delta+\mathcal{P}_{1}-\lambda\mathbb{I})=det(\Delta-\lambda\mathbb{I})det(\mathbb{I}+(\Delta-\lambda\mathbb{I})^{-1}\mathcal{P}_{1})=0,$
(27)

where $\Delta=\begin{bmatrix}A_{1}&B_{1}\\\ 0&B_{2}\end{bmatrix}$ is an upper
block triangular matrix with explicitly known eigenvalues, and
$\mathcal{P}_{1}=\begin{bmatrix}0&0\\\ A_{2}&0\end{bmatrix}$ is the rank-1
‘perturbation’. Note that this is the explicit form of the characteristic Eq.
25.

Next, we recall some results on Hamiltonian matrices and Riccati theory Abou-
Kandil  _et al._ (2012). A $2m\times 2m$ complex matrix $H$ is called
Hamiltonian if $JH+H^{*}J=0$, where $J=\begin{bmatrix}0&\mathbb{I}\\\
-\mathbb{I}&0\end{bmatrix}$. If $H$ is Hamiltonian, and $\lambda$ is an
eigenvalue of $H$, then $-\bar{\lambda}$ is also an eigenvalue of $H$.
Necessarily, every $2m\times 2m$ Hamiltonian matrix is of the form

$\displaystyle H=\begin{bmatrix}A&B\\\ C&D\end{bmatrix},$ (28)

where $D=-A^{*}$, and $B$ and $C$ are both Hermitian. Associated with $H$ is a
continuous time algebraic (matrix) Riccati equation (CARE): $XA-DX+XBX-C=0$.
We say that an $m\times m$ matrix $X$ that solves the CARE is a _stabilizing
solution_ if $A+BX$ is Hurwitz. We will need the following result:

######  Theorem 1.

(Chen et al.Chen and Huang (2018)) Let $H$ be a Hamiltonian matrix of Eq. 28,
with $A$ Hurwitz (i.e, all its eigenvalues of $A$ lie in the left half plane).
Suppose $H$ does not have any eigenvalues on the imaginary axis. Then, there
exists an orthogonal (‘Schur’) transformation $V$ such that

$\displaystyle V^{*}HV=\begin{bmatrix}H_{11}&H_{12}\\\ 0&H_{22}\end{bmatrix},$
(29)

where all eigenvalues of $H_{11}$ are the stable eigenvalues of $H$. If we
block partition $V=\begin{bmatrix}V_{11}&V_{12}\\\
V_{21}&V_{22}\end{bmatrix}$, then $\begin{bmatrix}V_{11}\\\
V_{21}\end{bmatrix}$ are the $m$ stable Schur vectors corresponding to the
stable block $H_{11}$. Furthermore, the matrix $V_{11}$ is invertible, and
$X_{+}=V_{21}V_{11}^{-1}$ is the unique Hermitian stablizing solution of CARE.

Note that the matrix $N$ in Eq. 26 is not Hamiltonian since $A_{2}$ is not
Hermitian. We consider its ‘symmetrized’ version,
$N_{s}=\begin{bmatrix}A_{1}&B_{1}\\\
\dfrac{A_{2}+A^{*}_{2}}{2}&B_{2}\end{bmatrix}$, which is a Hamiltonian matrix.
Using the definitions of the submatrices discussed above, it follows that the
characteristic equation of $N_{s}$ reduces to that of $N$, i.e., Eq. 27. Hence
$N$ and $N_{s}$ are similar matrices, i.e., they have the same eigenvalues.
This implies that there exists an invertible $2m\times 2m$ complex matrix $P$
s.t. $N_{s}=P^{-1}NP$. Note $P$ can be numerically obtained by eigenvalue
decomposition of $N$ and $N_{s}$. Then we have the following result:

######  Lemma 1.

Suppose $N$ defined in Eq. 26 does not have any eigenvalues on the imaginary
axis. Then, the BVP defined by Eq. 26 along with arbitrary initial condition
$Y_{1}(0)=Y_{10}$, and final condition
$\displaystyle\lim_{t\rightarrow\infty}Y_{2}(t)=0$ has a solution
$(Y_{1}(t),Y_{2}(t))$ s.t. $\displaystyle\lim_{t\rightarrow\infty}Y_{1}(t)=0$.

###### Proof.

Since $N_{s}$ and $N$ are similar, $N_{s}$ also doesn’t have any eigenvalues
on the imaginary axis. Then $N_{s}$ satisfies the conditions of Theorem 1,
since $A_{1}$ is Hurwitz. Hence there exists a Hermitian stabilizing solution
$X_{+}$ for CARE corresponding to $N_{s}$. Define
$U=\begin{bmatrix}\mathbb{I}&0\\\ X_{+}&\mathbb{I}\end{bmatrix}$, such that
$U^{-1}=\begin{bmatrix}\mathbb{I}&0\\\ -X_{+}&\mathbb{I}\end{bmatrix}$. Then,
$U^{-1}N_{s}U=\begin{bmatrix}A_{c}&B\\\ 0&-A_{c}^{*}\end{bmatrix}$, where
$A_{c}=A_{1}+B_{1}X_{+}$ is a stable matrix by Theorem 1, and we have used the
fact the $X_{+}$ solves the CARE. Now define new states $(Z_{1},Z_{2})$ via
$PU\begin{bmatrix}Z_{1}\\\ Z_{2}\end{bmatrix}=\begin{bmatrix}Y_{1}\\\
Y_{2}\end{bmatrix}$. Note that $PU$ is invertible since it is the product of
two invertible matrices. The BVP system of Eq. 26 in the new variables is

$\displaystyle\begin{bmatrix}\dot{Z}_{1}\\\
\dot{Z}_{2}\end{bmatrix}=U^{-1}P^{-1}NPU\begin{bmatrix}Z_{1}\\\
Z_{2}\end{bmatrix}=\begin{bmatrix}A_{c}&B_{1}\\\
0&-A_{c}^{*}\end{bmatrix}\begin{bmatrix}Z_{1}\\\ Z_{2}\end{bmatrix}.$ (30)

From the second component of Eqs. 30, we get
$Z_{2}(t)=e^{-A_{c}^{*}t}Z_{2}(0)$. Since $A_{c}$ is Hurwitz, $Z_{2}$ blows up
unless we pick $Z_{2}(0)=0$. This implies $Z_{2}(t)=0$ for all $t\geq 0$.
Using this in the first component of Eqs. 30, we obtain
$Z_{1}(t)=e^{A_{c}t}Z_{1}(0)$, and hence,
$\lim_{t\rightarrow\infty}Z_{1}(t)=0$. By invertibility of $PU$, the above two
results imply that $\displaystyle\lim_{t\rightarrow\infty}Y_{1}(t)=0$, and
$\displaystyle\lim_{t\rightarrow\infty}Y_{2}(t)=0$. ∎

Note that if $P=\begin{bmatrix}P_{11}&P_{12}\\\ P_{21}&P_{22}\end{bmatrix}$,
the initial condition $Z_{1}(0)$ above is defined by
$(P_{11}+P_{12}X_{+})Z_{1}(0)=Y_{10}$. To obtain uniqueness, we need the
additional assumption that $(P_{11}+P_{12}X_{+})$ is invertible, in which case
the unique solution is
$(Y_{1}(t)=(P_{11}+P_{12}X_{+})e^{A_{c}t}(P_{11}+P_{12}X_{+})^{-1}Y_{10},Y_{2}(t)=(P_{21}+P_{22}X_{+})e^{A_{c}t}(P_{11}+P_{12}X_{+})^{-1}Y_{10})$.

From the above results, it follows that if the linearized MFG system of Eq. 25
has no eigenvalues on the imaginary axis, and the invertibility assumption
holds, the equilibrium solution $(\rho_{\xi},h_{\xi})$ of the MFG system
(17,18) is linearly stable.

###  III.3 Numerical results

To illustrate our theoretical results, we fix parameters $L=10,h=5,\sigma=2$,
and compute the MFG eigenvalues (i.e., eigenvalues of matrix $N$) by
numerically solving the algebraic Eq. 25 using 22 eigenfunctions for each
Fourier mode, i.e., $0\leq p,q\leq 21$. As in the forward equation case, the
$k=1$ Fourier mode is the relevant spatial mode for studying stability. Fig. 3
shows the eigenspectrum for $k=1$ for various values of the unit control cost
$r$. Since the matrix $N$ has the same spectrum as the Hamiltonian matrix
$N_{s}$, this spectrum is symmetric about the imaginary axis. At $r=1.4$,
there are no imaginary axis eigenvalues. As $r$ is decreased from 1.4, a pair
of eigenvalues approaches the imaginary axis, and eventually collides on it at
$r_{c}\approx 0.95$. As $r$ is reduced futher, the two eigenvalue move away
from each other up/down the imaginary axis.

(a) r=1.4

(b) r=1.1

(c) r=0.95

(d) r=0.8

Figure 3: The spectrum (close to the imaginary axis) of the linearized MFG
operator in Eq. 23 as the control cost $r$ is varied, for Fourier mode $k=1$.

According to Lemma 1, this implies that for $r>r_{c}$, the MFG equilibrium
$(\rho_{\xi},h_{\xi})$ is linearly stable in the sense of Definition 1. We use
Schur decomposition to compute $X_{+}$ using the formula in Theorem 1, and
follow the construction of the BVP solutions $(Y_{1}(t),Y_{2}(t))$ in Lemma 1.
Figure 4 shows that the norm of the unique solution for an arbitrarily chosen
$Y_{10}$ for $r=1.4$ decays to zero. The critical value of unit control cost
decreases upon increasing the noise intensity $\sigma$, as shown in Fig. 5.

Figure 4: The norm of $Y_{1}(t)$ (bold) and $Y_{2}(t)$ (dashed) as a function
of time for the unique solution $(Y_{1}(t),Y_{2}(t))$ of the BVP of Eq. 26.
Here, $r=1.4>r_{c}$, and $k=1$. We choose an arbitrary $Y_{1}(0)$, and the
corresponding value of $Y_{2}(0)$ is assigned according to Lemma 1. Figure 5:
The critical unit control cost $r_{c}$ as a function of $\sigma$, for $h=5$.

#### Travelling wave solutions of the MFG

We look for solutions bifurcating from the spatially homogeneous equilibrium
$(\rho_{\xi},h_{\xi})$ for $r<r_{c}$. It is known that MFGs exhibit the
turnpike property Zaslavski (2005); Cardaliaguet  _et al._ (2013), i.e., the
solution over a large but finite time horizon spends most of its time near the
solution of the infinite time problem (the so-called ‘ergodic’ solution). We
implment a Picard-iteration based algorithm for MFGs Lauriere (2021) over a
large time-horizon in Dedalus, and find that for each $r$ in an open interval
$r\in(r_{c}-\epsilon,r_{c})$, the ergodic solution is a travelling wave
solution of the form $(\rho(x-\omega t,u),h(x-\omega t,u))$ to the MFG Eqs.
(17,18). Fig. 6 shows the marginal density of the travelling wave solution for
$r=0.8$. No travelling wave solutions were found for $r>r_{c}$, in which case
the algorithm always converged to one of the two stable equilibria
$(\rho_{\pm\xi^{*}},h_{\pm\xi^{*}})$.

Figure 6: The marginal density $\bigintssss_{\mathbb{R}}\rho(t,x,v)dv$ of a
travelling wave solution of the MFG Eqs. (17,18) for $r=0.8<r_{c}$.

##  IV Conclusions and Discussion

Via the MFG framework, we have studied transitions between distinct collective
behaviors in a population of non-cooperative motile inertial agents that are
interacting with their neighbors over a finite distance, and minimizing a
biologically inspired cost function. The kinetic MFG model is shown to mimic
phase transitions previously observed in the phenomenological Czirók model.

The linear stability of the equilibrium states is equivalent to the existence
of a unique decaying solution to the linearized BVP (in time) derived from the
nonlinear MFG PDE system. We provide conditions on the spectrum of the linear
operator for such a solution to exist. The explicit calculations are carried
out using Fourier-Hermite discretization of the linearized PDE, and use
properties of Hamiltonian matrices and Riccati equations. The existence of a
traveling wave solution of the MFG when the equilibrium loses stability is
shown numerically. A rigorous bifurcation analysis will be taken up in a
future work.

While non-equilibrium systems cannot be described by a variational principle,
the MFG inverse modelling approach adopted here is based on a generalized
optimality principle, which can potentially be extended to kinetic and
hydrodynamic descriptions of other systems with decision-making agents. We
plan to study some of these extensions in the near future.

## Acknowledgments

This material is based upon work supported by the US National Science
Foundation under Grant No. 2102112.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Physical review letters 75, 1226 (1995).
  * Levine, Rappel, and Cohen (2000) H. Levine, W.-J. Rappel, and I. Cohen, “Self-organization in systems of self-propelled particles,” Physical Review E 63, 017101 (2000).
  * Strogatz (2000) S. H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” Physica D: Nonlinear Phenomena 143, 1–20 (2000).
  * Cucker and Smale (2007) F. Cucker and S. Smale, “Emergent behavior in flocks,” IEEE Transactions on automatic control 52, 852–862 (2007).
  * Gompper  _et al._ (2020) G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, _et al._ , “The 2020 motile active matter roadmap,” Journal of Physics: Condensed Matter 32, 193001 (2020).
  * Kassabov, Strogatz, and Townsend (2021) M. Kassabov, S. H. Strogatz, and A. Townsend, “Sufficiently dense Kuramoto networks are globally synchronizing,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021).
  * Ding  _et al._ (2022) L. Ding, W. Li, S. Osher, and W. Yin, “A mean field game inverse problem,” Journal of Scientific Computing 92, 7 (2022).
  * Caines (2021) P. E. Caines, “Mean field games,” in  _Encyclopedia of systems and control_ (Springer, 2021) pp. 1197–1202.
  * Bensoussan, Frehse, and Yam (2013) A. Bensoussan, J. Frehse, and P. Yam,  _Mean field games and mean field type control theory_ , Vol. 101 (Springer, 2013).
  * Huang, Caines, and Malhame (2007) M. Huang, P. E. Caines, and R. P. Malhame, “Large-population cost-coupled LQG problems with nonuniform agents: Individual-mass behavior and decentralized Nash equilibria,” IEEE Transactions on Automatic Control 52, 1560–1571 (2007).
  * Solon  _et al._ (2015) A. P. Solon, J.-B. Caussin, D. Bartolo, H. Chaté, and J. Tailleur, “Pattern formation in flocking models: A hydrodynamic description,” Physical Review E 92, 062111 (2015).
  * Orosz, Krauskopf, and Wilson (2005) G. Orosz, B. Krauskopf, and R. E. Wilson, “Bifurcations and multiple traffic jams in a car-following model with reaction-time delay,” Physica D: Nonlinear Phenomena 211, 277–293 (2005).
  * Bardi and Fischer (2019) M. Bardi and M. Fischer, “On non-uniqueness and uniqueness of solutions in finite-horizon mean field games,” ESAIM: Control, Optimisation and Calculus of Variations 25, 44 (2019).
  * Nourian, Caines, and Malhamé (2011) M. Nourian, P. E. Caines, and R. P. Malhamé, “Mean field analysis of controlled Cucker-Smale type flocking: Linear analysis and perturbation equations,” IFAC Proceedings Volumes 44, 4471–4476 (2011).
  * Grover, Bakshi, and Theodorou (2018) P. Grover, K. Bakshi, and E. A. Theodorou, “A mean-field game model for homogeneous flocking,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 061103 (2018).
  * Cirant (2019) M. Cirant, “On the existence of oscillating solutions in non-monotone mean-field games,” Journal of Differential Equations 266, 8067–8093 (2019).
  * Cirant, Nurbekyan  _et al._ (2018) M. Cirant, L. Nurbekyan, _et al._ , “The variational structure and time-periodic solutions for mean-field games systems,” Minimax Theory and its applications 3, 227–260 (2018).
  * Ullmo, Swiecicki, and Gobron (2019) D. Ullmo, I. Swiecicki, and T. Gobron, “Quadratic mean field games,” Physics Reports 799, 1–35 (2019).
  * Salhab, Malhame, and Le Ny (2019) R. Salhab, R. P. Malhame, and J. Le Ny, “Collective stochastic discrete choice problems: A Min-LQG dynamic game formulation,” IEEE Transactions on Automatic Control (2019).
  * Hongler (2020) M.-O. Hongler, “Mean-field games and swarms dynamics in Gaussian and non-Gaussian environments.” Journal of Dynamics & Games 7 (2020).
  * Lori and Grover (2024) A. A. R. Lori and P. Grover, “Topological bifurcations in a mean-field game,” arXiv preprint arXiv:2405.05473 (2024).
  * Yin  _et al._ (2012) H. Yin, P. G. Mehta, S. P. Meyn, and U. V. Shanbhag, “Synchronization of coupled oscillators is a game,” IEEE Transactions on Automatic Control 57, 920–935 (2012).
  * Carmona, Cormier, and Soner (2023) R. Carmona, Q. Cormier, and H. M. Soner, “Synchronization in a Kuramoto mean field game,” Communications in Partial Differential Equations 48, 1214–1244 (2023).
  * Cesaroni and Cirant (2024) A. Cesaroni and M. Cirant, “Stationary equilibria and their stability in a Kuramoto MFG with strong interaction,” Communications in Partial Differential Equations , 1–27 (2024).
  * Kuramoto (1975) Y. Kuramoto, “International symposium on mathematical problems in theoretical physics,” Lecture notes in Physics 30, 420 (1975).
  * Czirók, Barabási, and Vicsek (1999) A. Czirók, A.-L. Barabási, and T. Vicsek, “Collective motion of self-propelled particles: Kinetic phase transition in one dimension,” Physical Review Letters 82, 209 (1999).
  * Degond, Liu, and Ringhofer (2014) P. Degond, J.-G. Liu, and C. Ringhofer, “Large-scale dynamics of mean-field games driven by local Nash equilibria,” Journal of Nonlinear Science 24, 93–115 (2014).
  * Barre and Métivier (2016) J. Barre and D. Métivier, “Bifurcations and singularities for coupled oscillators with inertia and frustration,” Physical Review Letters 117, 214102 (2016).
  * Acebrón, Bonilla, and Spigler (2000) J. Acebrón, L. Bonilla, and R. Spigler, “Synchronization in populations of globally coupled oscillators with inertial effects,” Physical Review E 62, 3437 (2000).
  * Laub (1991) A. J. Laub, “Invariant subspace methods for the numerical solution of Riccati equations,” in  _The Riccati Equation_ (Springer, 1991) pp. 163–196.
  * Abou-Kandil  _et al._ (2012) H. Abou-Kandil, G. Freiling, V. Ionescu, and G. Jank,  _Matrix Riccati equations in control and systems theory_ (Birkhäuser, 2012).
  * Huang and Zhou (2019) M. Huang and M. Zhou, “Linear quadratic mean field games: Asymptotic solvability and relation to the fixed point approach,” IEEE Transactions on Automatic Control 65, 1397–1412 (2019).
  * Garnier, Papanicolaou, and Yang (2019) J. Garnier, G. Papanicolaou, and T.-W. Yang, “Mean field model for collective motion bistability.” Discrete & Continuous Dynamical Systems-Series B 24 (2019).
  * Abramowitz and Stegun (1964) M. Abramowitz and I. Stegun,  _Handbook of Mathematical Functions with Formulas, Graph and Mathematical Tables_ (Dover, 1964).
  * Burns  _et al._ (2020) K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, “Dedalus: A flexible framework for numerical simulations with spectral methods,” Phys. Rev. Res. 2, 023068 (2020).
  * Nourian (2013) M. Nourian,  _Mean field game theory: consensus, leader-follower and major-minor agent systems_ , Ph.D. thesis, McGill University (2013).
  * Huang, Malhamé, and Caines (2006) M. Huang, R. P. Malhamé, and P. E. Caines, “Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle,” Communications in Information & Systems 6, 221–252 (2006).
  * Guéant (2009) O. Guéant, “A reference case for mean field games models,” Journal de mathématiques pures et appliquées 92, 276–294 (2009).
  * Chen and Huang (2018) X. Chen and M. Huang, “Linear-quadratic mean field control: the Hamiltonian matrix and invariant subspace method,” in  _2018 IEEE Conference on Decision and Control (CDC)_ (IEEE, 2018) pp. 4117–4122.
  * Zaslavski (2005) A. Zaslavski,  _Turnpike properties in the calculus of variations and optimal control_ , Vol. 80 (Springer Science & Business Media, 2005).
  * Cardaliaguet  _et al._ (2013) P. Cardaliaguet, J.-M. Lasry, P.-L. Lions, and A. Porretta, “Long time average of mean field games with a nonlocal coupling,” SIAM Journal on Control and Optimization 51, 3558–3591 (2013).
  * Lauriere (2021) M. Lauriere, “Numerical methods for mean field games and mean field type control,” Mean field games 78, 221 (2021).