# Stability of stationary states for mean field models with multichromatic
interaction potentials

Benedetta Bertoli1, Benjamin D. Goddard2, Grigorios A. Pavliotis1

###### Abstract

We consider weakly interacting diffusions on the torus, for multichromatic
interaction potentials. We consider interaction potentials that are not
$H$-stable, leading to phase transitions in the mean field limit. We show that
the mean field dynamics can exhibit multipeak stationary states, where the
number of peaks is related to the number of nonzero Fourier modes of the
interaction. We also consider the effect of a confining potential on the
structure of non-uniform steady states. We approach the problem by means of
analysis, perturbation theory and numerical simulations for the interacting
particle systems and the PDEs.

1Department of Mathematics, Imperial College London, London SW7 2AZ, UK.
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
2School of Mathematics and the Maxwell Institute for Mathematical Sciences,
University of Edinburgh, Edinburgh EH9 3FD, UK<EMAIL_ADDRESS>
##  1 Introduction

Nonlinear and nonlocal Fokker-Planck (advection-diffusion) equations appear in
several applications, including stellar dynamics [6], plasma physics [2],
mathematical biology [30, 33], active matter [32], biophysics [18][Sec 5.3],
nematic liquid crystals [11] and models for opinion formation [22]. Such PDEs
can exhibit nontrivial, i.e. non-uniform, stationary states, describing
collective behaviour and the emergence of coherent structures, as an effect of
interactions between agents at the microscale. In recent years, great progress
has been made towards the understanding of the emergence of such collective
behaviour. The purpose of this paper is to study the creation and stability of
multimodal/multipeak stationary states for nonlinear, nonlocal Fokker-Planck
equations on the torus.

In this paper, we will consider nonlinear, nonlocal Fokker-Planck equations of
the form

$\frac{\partial\rho}{\partial t}=\beta^{-1}\Delta\rho+\nabla\cdot(\nabla
V\rho)+\nabla\cdot((\nabla W\star\rho)\rho),\quad\rho(x,0)=\rho_{0}(x),$ (1)

on $\mathbb{T}^{d}:=[0,2\pi]^{d}$ with periodic boundary conditions. We will
focus on the one dimensional problem. Here $\rho$ denotes the
density/distribution function, $\rho_{0}$ the initial condition, $\beta$ the
inverse temperature and $V$ and $W$ the confining and (symmetric) interaction
potentials, respectively. Examples of dynamics described by (1) are the Haken-
Kelso-Bunz model from biophysics [18] with
$V(x)=-\alpha\cos(x)-\gamma\cos(2x)$ and $W(x)=-K\cos(x)$, where
$\alpha,\,\gamma,K$ are constants. We also mention the XY ($O(2)$) model with
an external magnetic field that was studied in [14], corresponding to
$V(x)=-\alpha\cos(x)$ and $W(x)=-K\cos(x)$, as well as the the noisy
Kuramoto/Brownian mean field model [5, 10], $V\equiv 0$ and $W(x)=-\cos(x)$
and the noisy Hegselmann-Krause model for opinion dynamics [22, 19]. Several
additional examples can be found in [9][Sec. 6]. A very nice presentation of
the nonlinear Fokker-Planck equation on the torus from a theoretical physics
perspective can be found in [18, Sec. 5.3] and in [10].

As is well known [9], the McKean-Vlasov PDE has a gradient structure:

$\partial_{t}\rho=\nabla\cdot\left(\rho\nabla\frac{\delta\mathcal{F}}{\delta\rho}\right),$
(2)

where $\mathcal{F}$ denotes the free energy

$\displaystyle\mathcal{F}(\rho)=\beta^{-1}\int_{\mathbb{T}^{d}}\rho(x)\log(\rho(x))\,\mathrm{d}x+\int_{\mathbb{T}^{d}}V(x)\rho(x)\,\mathrm{d}x+\frac{1}{2}\int_{\mathbb{T}^{d}}\int_{\mathbb{T}^{d}}W(x-y)\rho(x)\rho(y)\,\mathrm{d}x\,\mathrm{d}y.$
(3)

Stationary states of the mean field dynamics can be characterized as critical
points of the free energy [9]. The main goal of this paper is to study the
dynamical stability of such states, in particular of stationary states that
describe collective, organized behaviour.

Stationary states of the McKean-Vlasov PDE satisfy the Kirkwood-
Monroe/generalized Lane-Emden integral equation [4, 9, 18, 25]

$\rho=\frac{1}{Z}e^{-\beta(V+W\star\rho)},\quad
Z=\int_{\mathbb{T}^{d}}e^{-\beta(V+W\star\rho)}\,\,\mathrm{d}x.$ (4)

In the absence of an external potential, the uniform distribution, describing
the disordered state, is always a stationary state of the Fokker-Planck
equation (1). Collective, organized behaviour, described by localized or
multipeak solutions, becomes possible when the disordered state becomes
unstable.

As is well-known, the Fokker-Planck equation (1), arises in the mean field
limit of a system of weakly interacting diffusions [8, 29, 10, 26]. In
particular, we consider a system of interacting diffusions of the form:

$\,\mathrm{d}x_{i}(t)=-\nabla
V(x_{i})\,\mathrm{d}t-\frac{1}{N}\sum_{j=1}^{N}\nabla
W(x_{i}-x_{j})\,\mathrm{d}t+\sqrt{2\beta^{-1}}\,\mathrm{d}B_{i}(t),\quad
x_{i}(0)\sim\rho_{0},$ (5)

for $1\leq i\leq N$, where $B_{i}(t)$ denote standard $d$-dimensional
independent Brownian motions. Under appropriate assumptions on the confining
and interaction potentials, and for chaotic initial conditions, the sequence
of empirical measures $\rho^{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{i}(t)}$
converges to the solution of the mean field PDE (1). Equivalently, the
$N$-particle distribution function for the interacting particle system (5) can
be written as $\rho^{N}(x_{1},\dots
x_{N},t)\approx\Pi_{j=1}^{N}\rho(x_{j},t)$. Rigorous convergence results,
either at the level of the empirical measure or of the product measure
structure of the $N$-particle distribution function, are by now well-
estabilished and we refer to, e.g. [29, 24, 8]. We will refer to this equation
as the McKean-Vlasov PDE.

###  1.1 Literature Review

There is extensive literature on the calculation of stationary states for the
McKean-Vlasov PDE and on the study of their stability as well as on
applications to mathematical biology, in particular mass-selection in
alignment models with non-deterministic effects and in active matter. The
number of stationary states of the McKean-Vlasov PDE and their stability has
been studied extensively, either by studying the Kirkwood-Monroe map [4], or
by studying critical points of the free energy functional [25] or by studying
the stationary McKean-Vlasov PDE.111All these approaches are, of course,
equivalent; see [9, Prop. 2.4]. The existence and stability of multipeak
solutions, the problem that we will primarily focus on in this paper, was
studied in [33, 20] using PDEs/ODEs techniques. In [33], the existence of a
2-peak steady state is proved under appropriate assumptions on the interaction
potential. The stability of steady states was investigated numerically; the
simulations presented in this paper suggest that only $1$ and $2$-peak steady
states can be stable, while solutions with $4$ peaks are always unstable. The
work in [20] builds on these results. Bifurcation theory for the stationary
McKean-Vlasov equation for multichromatic interaction potentials was recently
analysed in [35]; in this paper the Hodgkin-Huxley oscillator model is
studied, with an interaction potential consisting of two Fourier modes with
opposite sign, $W(\theta)=-\cos(\theta)+\varepsilon\cos(2(\theta-\alpha))$,
$\varepsilon\geq 0,\alpha\in(-\pi/2,\pi/2)$. Critical points of the free
energy functional for the Onsager model for liquid crystals, $W(x)=|\sin(x)|$
have been studied extensively. See, e.g. [25, 17, 16] and the references
therein.

A quite comprehensive theory of bifurcations from the uniform distribution and
of phase transitions for the McKean-Vlasov PDE on the torus was developed in
[9]. The goal of the present paper is to study, by means of analysis,
systematic perturbation theory and numerical simulations, the stability of
non-uniform states, and in particular of multipeak solutions. The stability of
the non-uniform state for the noisy Kuramoto model was studied using spectral
theoretic arguments in [5]. One of the goals of the present study is to extend
the analysis from this paper to multichromatic interaction potentials.

###  1.2 Our Contributions

In this paper we consider the stability of multipeak stationary solutions for
the one-dimensional McKean-Vlasov PDE, both in the presence or absence of a
confining potential. Our main contributions are the following:

  * •

We provide a detailed stability analysis of the non-uniform state for
multichromatic potentials.

  * •

In particular, we find that the uniform state changes stability at some
critical value of the temperature $\beta$, and that multi-peak stationary
states are unstable.

  * •

We calculate the eigenvalues of the linearized McKean-Vlasov operator above
the bifurcation.

  * •

We present very detailed numerical experiments by solving the evolution PDE,
the SDEs for the interacting particle system, and the eigenvalue problem for
the linearized operator.

The rest of the paper is organized as follows. In Section 2 we present the
models that we will consider in this paper and we analyse the self-consistency
equation(s). In Section 3 we study the stability of stationary states by
either calculating the second variation of the free energy or by linearizing
the McKean-Vlasov PDE. Peturbative results for the eigenvalues of the
linearized McKean-Vlasov operator close to the bifurcation point are presented
in Section 4. The results of extensive numerical simulations based on both the
PDE and SDE formulations are shown in Section 5. Conclusions and comments on
future work are presented in Section 6.

##  2 Set-up and self-consistency equations

###  2.1 Phase transitions, stability analysis and the self-consistency
equation

For $H$-stable potentials, i.e. interaction potentials with non-negative
Fourier coefficients, and in the absence of a confining potential, the free
energy functional is convex and the uniform distribution is the unique,
globally stable stationary state [4].

For monochromatic interaction potentials of the form
$W(x)=-\cos(kx),k\in\mathbb{N}$, and in the absence of a confining potential,
a detailed characterization of stationary states was given in [9]. In
particular, we have the following:

######  Proposition 2.1.

The generalised Kuramoto model $W(x)=-\cos(kx)$, for some $k\in\mathbb{N}$,
$k\neq 0$ exhibits a continuous transition point at the linear instability
threshold $\kappa_{c}=\kappa_{\sharp}$. Additionally, for $\kappa>\kappa_{c}$,
the equation $F(\varrho,\kappa)=0$ has only two solutions in $L^{2}(U)$ (up to
translations). The nontrivial one, $\varrho_{\kappa}$ minimises
$\mathcal{F}_{\kappa}$ for $\kappa>\kappa_{c}$ and converges in the narrow
topology as $\kappa\to\infty$ to a normalised linear sum of equally weighted
Dirac measures centred at the minima of $W(x)$.

We will consider even confining and interaction potentials with a finite
number of non-zero Fourier modes:

$\displaystyle V(x)=\sum_{k=1}^{m}v_{k}\cos(kx),\quad\mbox{and}\quad
W(x)=\sum_{k=1}^{n}a_{k}\cos(kx),$

with $\\{v_{k}\\}_{k=1}^{m}$ real-valued and $\\{a_{k}\\}_{k=1}^{n}$ non-
positive.

The main objective of this paper is to study the nature and stability of
stationary states of the McKean-Vlasov PDE by analysing the Kirkwood-Monroe
integral equation (4). We reiterate that, for nontrivial confining potentials
$V(x)$, the uniform state is no longer a stationary state. In particular, if
$W(x)=\sum_{k=1}^{n}a_{k}\cos(kx)$, $V(x)=\sum_{k=1}^{m}v_{k}\cos(kx)$, then
the stationary state is:

$\displaystyle\rho=\frac{1}{Z}\exp\left(\beta\sum_{k=1}^{\max\\{m,n\\}}(v_{k}+r_{k}a_{k})\cos(kx)\right),$

where
$Z=\int_{0}^{2\pi}\exp\left(\beta\sum_{k=1}^{\max\\{m,n\\}}(v_{k}+r_{k}a_{k})\cos(kx)\right)$,
and we set $v_{k}=0$ for $k>m$, and $a_{k}=0$ for $k>n$.

The Fourier coefficients of $\rho$,
$r_{k}=\int_{0}^{2\pi}\cos(kx)\rho(x)\,\mathrm{d}x$, solve the self-
consistency equations:

$\displaystyle
r_{l}=\frac{1}{Z}\int_{0}^{2\pi}\cos(lx)\exp\left(\beta\sum_{k=1}^{\max\\{m,n\\}}(v_{k}+r_{k}a_{k})\cos(kx)\right),$

for $1\leq l\leq\max\\{m,n\\}$

These self-consistency equations, in the absence of a confining potential,
were first derived in [3].

We are interested in the number of critical points of the steady states of the
McKean-Vlasov equation (1). When there is no external potential, it is
conceivable that one may relate this to the number of Fourier modes of the
interaction potential. However, as we show below, even if one fixes $W$, then
this is not possible for general external potentials. As a particular example,
one may fix $W$ and (smooth, non-negative) $\rho$ and then choose $V$ such
that $\rho$ is a steady state, irrespective of the number of peaks it
contains.

######  Lemma 2.2.

Let $\hat{\rho}$ be a smooth, non-negative, $2\pi$-periodic function and fix
the interaction kernel $W$. Then there exists an external potential, $V$, such
that $\hat{\rho}$ is an equilibrium of the McKean-Vlasov equation (1).

###### Proof.

By the non-negativity of $\hat{\rho}$, we may write
$\hat{\rho}(x)=\exp\big{(}-\beta f(x)\big{)}$ for some $f(x)$, which is fixed
by the choice of $\hat{\rho}$. Similarly, by (4), if
$\hat{\rho}(x)=\exp\big{[}-\beta\big{(}V+W\star\rho-\ln(Z)\big{)}\big{]}$,
then $\hat{\rho}$ is an equilibrium of (1). Note that, for convenience, we
have written the normalization constant in the exponent. Hence it is clear
that we want $f(x)=V(x)+(W\star\rho)(x)-\ln(Z)$, which can be achieved by
taking $V(x)=f(x)+(W\star\rho)(x)-c$, where $c$ is a constant chosen to ensure
the correct normalization. ∎

We note that, in fact, the maximum number of non-zero Fourier modes required
to produce an equilibrium distribution with $M$ non-zero Fourier modes under a
kernel with $n$ non-zero Fourier modes is $\max\\{n,M\\}$.

There are many well-known results for the interaction potential
$W(x)=-\cos(x)$, which gives rise to the Kuramoto model. In particular, it is
known (see [5]) that the Fokker-Planck PDE for the Kuramoto model undergoes a
phase transition. This means that there exists a critical value $\beta_{c}$ of
the inverse temperature $\beta$ such that, for $\beta<\beta_{c}$, the PDE
admits a unique stationary solution (the uniform state $1/(2\pi)$), while for
$\beta>\beta_{c}$, the PDE has multiple stationary states. Each of these
stationary solutions can be written as $q(x)=\rho(x+x_{0})$ for some
$x_{0}\in(0,2\pi)$, where

$\displaystyle\rho(x)=\frac{1}{Z}\exp(\beta r\cos(x)),$

with $Z=\int_{0}^{2\pi}\exp(\beta r\cos(y))\,\mathrm{d}y$. Here $r$ is a
solution of the self-consistency equation:

$\displaystyle r=\frac{\int_{0}^{2\pi}\cos(x)\exp(\beta
r\cos(x))\,\mathrm{d}x}{\int_{0}^{2\pi}\exp(\beta r\cos(x))\,\mathrm{d}x}.$
(6)

The number of stationary solutions is then given by the number of solutions
$r$ to the self-consistency equation. As an example in which the self-
consistency equation can be solved analytically, we consider the following
Brownian mean field model in a magnetic field studied in [14], for the
confining and interaction potentials $V(x)=-\eta\cos(x),W(x)=-a\cos(x)$,
example from Michela’s paper. This can be written as an unbounded spin system
with Hamiltonian, with ${\bf s}_{i}=(\cos(x),\sin(x))$:

$H=-\sum_{i=1}^{N}{\bf J}\cdot{\bf s}_{i}-\frac{2}{N}\sum_{i,j=1}^{N}{\bf
s}_{i}\cdot{\bf s}_{j},$ (7)

where ${\bf J}=-\eta(1,0)$. In this case we can calculate the stationary
states analytically [14]: for $\beta<\beta_{c}$ there exists a unique steady
state

$\rho_{min}(x)=\frac{1}{Z_{min}}e^{a^{min}\cos(2\pi x)},\quad
Z_{min}=\int_{\mathbb{T}}e^{a^{min}\cos(2\pi x)\,dx},$ (8)

Furthermore, for $\beta>\beta_{c}$ there exist at least two steady states; in
addition to (8), we have

$\rho_{*}(x)=\frac{1}{Z_{min}}e^{a^{*}\cos(2\pi x)},\quad
Z_{*}=\int_{\mathbb{T}}e^{a^{*}\cos(2\pi x)\,dx},$ (9)

for some $a^{*}=a^{*}(\beta)$, with $a^{*}<0<a^{min}$. $\rho^{min}$ is the
unique minimizer and $\rho^{*}$ is a non-minimizing critical point of the
periodic mean field energy.

Our aim is to extend this type of result to more general interaction
potentials $W$ with a varying number of Fourier modes, such as the bichromatic
potential $W(x)=a_{1}\cos(x)+a_{2}\cos(2x)$, $a_{1},a_{2}<0$. We address the
question of existence and uniqueness of stationary states and their stability.

Following the same calculations as the ones done in [5] for the Kuramoto
model, we find that invariant measures for the corresponding Fokker-Planck PDE
are given by:

$\rho(x)=\frac{\exp(\beta(a_{1}r_{1}\cos(x)+a_{2}r_{2}\cos(2x)))}{\int_{0}^{2\pi}\exp(\beta(a_{1}r_{1}\cos(y)+a_{2}r_{2}\cos(2y)))\,\mathrm{d}y},$
(10)

with $r_{1}$ and $r_{2}$ satisfying the self-consistency equations:

$\displaystyle
r_{1}=\frac{\int_{0}^{2\pi}\cos(x)\exp(\beta(a_{1}r_{1}\cos(x)+a_{2}r_{2}\cos(2x)))\,\mathrm{d}x}{\int_{0}^{2\pi}\exp(\beta(a_{1}r_{1}\cos(y)+a_{2}r_{2}\cos(2y)))\,\mathrm{d}y},$

and

$\displaystyle
r_{2}=\frac{\int_{0}^{2\pi}\cos(2x)\exp(\beta(a_{1}r_{1}\cos(x)+a_{2}r_{2}\cos(2x)))\,\mathrm{d}x}{\int_{0}^{2\pi}\exp(\beta(a_{1}r_{1}\cos(y)+a_{2}r_{2}\cos(2y)))\,\mathrm{d}y}.$

We will prove this result for the more general $n$-modes potential in the next
section.

In this case, it is harder to rigorously deduce results about $r_{1}$ and
$r_{2}$. However, as this potential only induces two self-consistency
equations, we can solve these numerically. We take, for example,
$a_{1}=a_{2}=-1$. The number of solutions of these equations will determine
the number of stationary distributions of the Fokker-Planck PDE. It is easy to
notice that for any value of $\beta$, $r_{1}=r_{2}=0$ is a solution, which
corresponds to the uniform stationary state $\rho=1/(2\pi)$. Numerical
experiments indicate that for $\beta<2$, this is the only solution. However,
when $\beta>2$, we start seeing two other pairs of solutions $(r_{1},r_{2})$;
substituting these into (10) gives us three different steady states. In Figure
1 we show these steady states for for $\beta=3$ and $\beta=10$, along with the
corresponding values for $r_{1}$ and $r_{2}$.

Figure 1: Probability distributions corresponding to the three different
solutions for $(r_{1},r_{2})$ for $\beta=3$ (left) and $\beta=10$ (right).

Note that, alongside the uniform distribution, there is a one-peak steady
state (as also seen for the Kuramoto model), and we now have a further
solution with two peaks, reflecting the multichromatic nature of the
interaction potential.

In the following section we aim to address the question of stability of such
solutions.

##  3 Stability Analysis

We perform the stability calculations in two ways. As before, for interaction
potentials that are not $H-$stable we expect a critical value of $\beta$ above
which the McKean-Vlasov PDE admits stationary solutions other than the uniform
state. This critical value is given by $\beta_{c}=\frac{1}{T_{c}}$, where
$T_{c}=\min_{1\leq k\leq n}\left\\{-\frac{a_{k}}{2}:1\leq k\leq n\right\\}$.
We will identify this critical temperature in three ways: firstly we analyse
the second variation of the free energy, then we perform a linear stability
analysis of the PDE, and finally we analyse the problem numerically. These
methods will also give us an insight on the stability of the states we find.

###  3.1 Second variation of the free energy

The first approach we use to perform the linear stability analysis and to
identify the continuous phase transition, based on the one used in [10],
consists in analysing the eigenvalues of the second variation of the free
energy. We first present the calculation in the absence of a confining
potential. We recall that the free energy of this system is given by

$\displaystyle
F(\rho)=\frac{1}{2}\int_{0}^{2\pi}\int_{0}^{2\pi}W(x-y)\rho(x)\rho(y)\,\mathrm{d}x\,\mathrm{d}y+\beta^{-1}\int_{0}^{2\pi}\rho(x)\log(\rho(x))\,\mathrm{d}x.$

Its second variation is

$\displaystyle\delta^{2}F[\delta\rho,\,\delta\rho]=\frac{1}{2}\int\delta\rho\delta\phi\,\mathrm{d}x+\frac{\beta^{-1}}{2}\int\frac{(\delta\rho)^{2}}{\rho}\,\mathrm{d}x,$

where $\delta\rho$ is the perturbation and:

$\displaystyle\delta\phi(x)=\int_{0}^{2\pi}W(x-y)\delta\rho(y)\,\mathrm{d}y.$

For our interaction potential $W(x)=-\sum_{k=1}^{n}a_{k}\cos(kx)$, this is
equal to:

$\displaystyle\delta^{2}F=(2\pi)^{2}\sum_{j=1}^{\infty}\left(\beta^{-1}+\frac{a_{j}}{2}\right)|\delta\hat{\rho}_{j}|^{2}.$

We then write:

$\displaystyle\delta\rho=\frac{\,\mathrm{d}q}{\,\mathrm{d}x},\text{ where
}q(x)=\int_{0}^{x}\delta\rho(y,t)\,\mathrm{d}y,$

to obtain:

$\displaystyle\delta^{2}F$
$\displaystyle=\frac{1}{2}\int_{0}^{2\pi}\delta\rho\delta\phi\,\mathrm{d}x+\frac{\beta^{-1}}{2}\int_{0}^{2\pi}\frac{(\delta\rho)^{2}}{\rho}\,\mathrm{d}x.$

This gives, for the first integral:

$\displaystyle\frac{1}{2}\int_{0}^{2\pi}\delta\rho\delta\phi\,\mathrm{d}x=\frac{1}{2}\int_{0}^{2\pi}\int_{0}^{2\pi}q(x)\frac{\,\mathrm{d}}{\,\mathrm{d}x}\frac{\,\mathrm{d}}{\,\mathrm{d}y}W(x-y)q(y)\,\mathrm{d}x\,\mathrm{d}y,$

and for the second integral:

$\displaystyle\frac{\beta^{-1}}{2}\int_{0}^{2\pi}\frac{(\delta\rho)^{2}}{\rho}\,\mathrm{d}x$
$\displaystyle=-\frac{\beta^{-1}}{2}\int_{0}^{2\pi}q(x)\frac{\,\mathrm{d}}{\,\mathrm{d}x}\left(\frac{1}{\rho(x)}\frac{\,\mathrm{d}}{\,\mathrm{d}x}q(x)\right)\,\mathrm{d}x.$

Therefore, the second variation of the free energy can be written as a
quadratic form/integral operator:

$\displaystyle\delta^{2}F[q,q]=\int_{0}^{2\pi}\int_{0}^{2\pi}K(x,y)q(x)q(y)\,\mathrm{d}x\,\mathrm{d}y,$

where the kernel is defined as

$K(x,y)=\frac{1}{2}\sum_{k=1}^{n}k^{2}a_{k}\cos(k(x-y))-\frac{\beta^{-1}}{2}\delta(x-y)\frac{\,\mathrm{d}}{\,\mathrm{d}x}\left(\frac{1}{\rho}\frac{\,\mathrm{d}}{\,\mathrm{d}x}\right).$

We are therefore led to considering the eigenvalue problem:

$2\pi\frac{\,\mathrm{d}^{2}q}{\,\mathrm{d}x^{2}}+\beta\int_{0}^{2\pi}q(y)\left(\sum_{k=1}^{n}k^{2}a_{k}\cos(k(x-y))\right)\,\mathrm{d}y=-2\lambda
q.$ (11)

The eigenmodes are $q_{m}=A_{m}\cos(mx)$, $q_{m}=B_{m}\sin(mx)$, where
$A_{m},B_{m}\in\mathbb{R}$. It is sufficient to consider the even
eigenfunctions. Using trigonometric identities, we can compute the eigenvalues
of the integral operator. We first conclude that, for $n<|m|$, we have
$\lambda_{m}=\pi m^{2}>0$, so these modes do not induce instability. For
$|m|\leq n$, we have

$\displaystyle\lambda_{m}=-\frac{1}{2}\left(-2\pi m^{2}-\beta
m^{2}a_{m}\pi\right)=\pi m^{2}\left(2+\beta a_{m}\right)$

This is positive for $\beta a_{m}<-2$, so $\beta<-\frac{a_{m}}{2}$. As we are
concerned with the first point of linear instability, the critical temperature
is $\beta_{c}:=\min_{k}\\{-\frac{2}{a_{k}}:1\leq k\leq n\\}$, see also [9]. As
an example, we note that by setting $a_{1}=-1,a_{k}=0$ for $2\leq n$, we
obtain the known result for the Kuramoto model, $\beta_{c}=2$. Similarly, for
the interaction potentials $W(x)=-3\cos(x)-\cos(2x)$ or
$W(x)=-\cos(x)-\cos(2x)-3\cos(3x)$, the critical temperature is
$\beta_{c}=\frac{2}{3}$, as the critical value is not influenced by which
Fourier mode has the highest coefficient, but rather by the magnitude of the
coefficient itself. Critical temperatures for more examples can be found,
together with the corresponding numerical simulations, in Section 5.

Our goal is now to extend this method to study the stability of nonuniform
steady states. This extension is not straightforward, and the details will be
presented in future work. Here we present only a specific example. We consider
$W(x)=-\cos(x)-\frac{1}{2}\cos(2x)$, $V(x)=0$. The nonuniform stationary state
is given by
$\tilde{\rho}=\frac{1}{Z}\exp\left(\beta\left(-r_{1}\cos(x)-\frac{r_{2}}{2}\cos(2x)\right)\right)$
with $r_{1},r_{2}\neq 0$. The eigenvalue problem for the second variation of
the free energy, computed at the nonuniform steady state, becomes:

$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}x}\left(\frac{1}{\tilde{\rho}}\frac{\,\mathrm{d}q}{\,\mathrm{d}x}\right)-\beta\int_{0}^{2\pi}q(y)(\cos(x-y)+2\cos(2(x-y))\,\mathrm{d}y=-2\lambda
q.$

We need to consider the following operator:

$\displaystyle\frac{\,\mathrm{d}^{2}q}{\,\mathrm{d}x^{2}}\left(Ze^{\beta
f(x)}\right)+\frac{\,\mathrm{d}q}{\,\mathrm{d}x}\left(Z\beta
f^{\prime}(x)e^{\beta
f(x)}\right)-\beta\int_{0}^{2\pi}q(y)\left(\cos(x-y)+4\cos(2(x-y))\right)\,\mathrm{d}y$

where $f(x):=-r_{1}\cos(x)-r_{2}\cos(2x)$.

We study the eigenvalue problem numerically for fixed values of $\beta$ and
$r_{1}$ and $r_{2}$. Let us first consider the pair $(r_{1},r_{2})$ that
appears immediately above the critical temperature, i.e. the pair for which
$r_{1},r_{2}\neq 0$, and which corresponds to the one-peak steady state. In
this case, we identify a change in the behaviour of the eigenvalues at around
$\beta=4.5$. This is the point at which the self consistency equations go from
having two solutions to three. For $\beta<4.5$, all eigenvalues are negative,
while at $\beta=4.5$ we start seeing some positive eigenvalues as well,
indicating that the stability of this steady state changes as the third
stationary distribution appears. On the other hand, the operator corresponding
to the multipeak solution (i.e. to the pair of coefficients $(0,r_{2})$) seems
to have a mix of positive and negative eigenvalues for any value of $\beta$.

###  3.2 Linearisation of the Fokker-Planck equation

We can calculate the value of the critical temperature also by looking at the
Fourier modes of the linearisation of the Fokker-Plank PDE:

$\displaystyle\frac{\partial\rho}{\partial t}(t,x)=\frac{\partial}{\partial
x}\left[\left(\int\rho(t,x-y)W^{\prime}(y)\,\mathrm{d}y\right)\rho(t,x)\right]+\beta^{-1}\frac{\partial^{2}\rho}{\partial
x^{2}}(t,x).$

Following the method in [19], we decompose
$\rho=\rho_{0}+\rho_{1}=\frac{1}{2\pi}+\rho_{1}$, where $\rho_{1}$ is a small
perturbation of $\rho$ so that $O(\rho_{1}^{2})$ is negligible. In the Fourier
domain the PDE becomes:

$\displaystyle\frac{\partial\hat{\rho}_{1}(t,j)}{\partial
t}=\left[i\rho_{0}j\int
e^{-ijy}W^{\prime}(y)\,\mathrm{d}y-\beta^{-1}j^{2}\right]\hat{\rho_{1}}(t,j),$

where $\hat{\rho}_{1}(j)$ is the $j$-th Fourier coefficient of $\rho_{1}$.
Therefore, the Fourier modes have growth rates:

$\displaystyle\gamma_{j}=\text{Re}\left[i\rho_{0}j\int
e^{-ijy}W^{\prime}(y)\,\mathrm{d}y-\beta^{-1}j^{2}\right]=\rho_{0}j\int\sin(jy)W^{\prime}(y)\,\mathrm{d}y-\beta^{-1}j^{2}.$

Substituting $W^{\prime}(x)=-\sum_{k=1}^{n}ka_{k}\sin(kx)$, our problem
reduces to studying the sign of:

$\displaystyle\gamma_{j}=-\frac{j}{2\pi}\int_{0}^{2\pi}\sin(jy)\sum_{k=1}^{n}ka_{k}\sin(ky)\,\mathrm{d}y-\beta^{-1}j^{2},$

as $\beta$ varies. For $j>n$, the integral term above is equal to $0$. Hence:

$\displaystyle\gamma_{j}=-\beta^{-1}j^{2}.$

As this is always negative, these modes do not induce instability. For $j\leq
n$:

$\displaystyle\gamma_{j}=-\frac{1}{2\pi}\int_{0}^{2\pi}\sin(jy)\sum_{k=1}^{n}ka_{k}\sin(ky)\,\mathrm{d}y-\beta^{-1}j=-\frac{ja_{j}}{2}-\beta^{-1}j.$

This is positive for $2\beta^{-1}<-a_{j}$, i.e. for $\beta>-\frac{2}{a_{j}}$,
as expected from the calculations in the previous sections.

##  4 Perturbation analysis

Our next goal is to study the stability of non-uniform stationary states close
to the critical interaction strength. One approach would be to linearize the
McKean-Vlasov operator

$\mathcal{L}\rho=\beta^{-1}\partial_{x}^{2}\rho+\partial_{x}\left(\partial_{x}V\rho\right)+\theta\partial_{x}\left((\partial_{x}W\star\rho)\rho\right)$
(12)

around the (non-uniform) stationary state $\rho_{\infty}$, the solution of the
stationary McKean-Vlasov PDE

$\beta^{-1}\partial_{x}^{2}\rho_{\infty}+\partial_{x}\left(\partial_{x}V\rho_{\infty}\right)+\theta\partial_{x}\left((\partial_{x}W\star\rho_{\infty})\rho_{\infty}\right)=0,$
(13)

on $[0,2\pi]$ with periodic boundary conditions, and where $\theta$ denotes
the interaction strength. The linearized operator is

$\displaystyle\mathcal{L}_{\rho_{\infty}}\rho$ $\displaystyle=$
$\displaystyle\beta^{-1}\partial_{x}^{2}\rho+\partial_{x}(\partial_{x}V\rho)+\theta\partial_{x}((\partial_{x}W\star\rho_{\infty})\rho)+\theta\partial_{x}((\partial_{x}W\star\rho)\rho_{\infty})$
(14) $\displaystyle=:$
$\displaystyle\beta^{-1}\partial_{x}^{2}\rho+\partial_{x}(\mathcal{U}_{\theta}\rho)+\int\kappa_{\theta}(x,y)\rho(y)\,dy.$

In writing the above, we have introduced the effective potential
$\mathcal{U}_{\theta}$. This is an integro-differential Fokker–Planck
operator. Via the standard ground state transformation we can map it to a
nonlocal Schrödinger operator [31][Sec. 4.9], of the form considered in [12,
13], albeit with different boundary conditions. The spectral properties of
(14) will be studied elsewhere.

Instead of considering (14), we follow [5][Sec. 2.5] (Eqn (2.57)) and perform
the ground state transformation _before_ the linearization to map the McKean-
Vlasov operator to a nonlinear and nonlocal Schrödinger operator. The
nonlinearity and nonlocality enter through the dependence on the Fourier modes
of the stationary states that satisfy the self-consistency equations. See also
[34, Sec. 3].

The Schrödinger operator for the noisy Kuramoto model is:

$\displaystyle
Hf=f^{\prime\prime}-\frac{1}{2}\left(\frac{\beta^{2}}{2}r^{2}\sin^{2}(x)+\beta
r\cos(x)\right)f.$

We consider the self-consistency equation (6) close to the critical inverse
temperature $\beta=2(1+\varepsilon$), for a small $\varepsilon>0$. For $\beta$
close to $2$, we use the approximation
$r(\beta)\approx\sqrt{1-\frac{2}{\beta}}$ [5][Sec. 2.5]. Substituting these
values of $r$ and $\beta$ in the expression for $H$:

$\displaystyle Hf$
$\displaystyle=f^{\prime\prime}-\left(\varepsilon\sin^{2}(x)+\varepsilon^{2}\sin^{2}(x)+\sqrt{\varepsilon(1+\varepsilon)}\cos(x)\right)f.$

Now, using the Taylor expansion for $\sqrt{1+\varepsilon}$, we have
$\sqrt{\varepsilon(1+\varepsilon)}=\sqrt{\varepsilon}+O(\varepsilon\sqrt{\varepsilon})$.
We set $\delta=\sqrt{\varepsilon}$ and can then rewrite the equation for $H$
as:

$\displaystyle
Hf=f^{\prime\prime}-\delta\cos(x)f-\delta^{2}\sin^{2}(x)f+O(\delta^{3})=:H_{0}+\delta
H_{1}+\delta^{2}H_{2}+O(\delta^{3}).$

Therefore, we can consider $H$ as a small perturbation of the operator
$H_{0}f=f^{\prime\prime}$. Our goal is to calculate the eigenvalues of $H$
perturbatively for small $\delta$. We consider the eigenvalue problem

$H\psi=E\psi.$ (15)

We expand $E$ and $\psi$ in power series in $\delta$:

$\displaystyle\psi=\psi_{0}+\delta\psi_{1}+\delta^{2}\psi_{2}+\ldots$
$\displaystyle E=E_{0}+\delta E_{1}+\delta^{2}E_{2}+\ldots$

We substitute these expansions into (15) to obtain the following sequence of
equations.

Order $O(1)$: $H_{0}\psi_{0}=E_{0}\psi_{0}$.

The eigenvalues of $H_{0}$ are given by $E^{(m)}=-m^{2}$ for $m\in\mathbb{N}$,
and the corresponding eigenfunctions are
$\psi^{(m)}(x)=A_{m}\cos(mx)+B_{m}\sin(mx)$, $A_{m},B_{m}\in\mathbb{R}$. Due
to the symmetry of the interaction potential, we will only consider even
eigenfunctions, i.e. $B_{m}=0$ for all $m\in\mathbb{N}$. We take
$A_{m}=\frac{1}{\sqrt{\pi}}$ so that $\langle\psi_{0},\psi_{0}\rangle=1$.

Order $O(\delta)$: $H_{0}\psi_{1}+H_{1}\psi_{0}=E_{0}\psi_{1}+E_{1}\psi_{0}$.

We take the inner product of both sides with $\psi_{0}$ to obtain:

$\displaystyle\langle\psi_{0},H_{0}\psi_{1}\rangle+\langle\psi_{0},H_{1}\psi_{0}\rangle=E_{0}\langle\psi_{0},\psi_{1}\rangle+E_{1}\langle\psi_{0},\psi_{0}\rangle.$

Since $H_{0}$ is self-adjoint in $L^{2}(0,2\pi)$, we obtain:

$\displaystyle
E_{1}=\frac{\langle\psi_{0},H_{1}\psi_{0}\rangle}{\langle\psi_{0},\psi_{0}\rangle}=0,$

Order $O(\delta^{2})$:
$H_{0}\psi_{2}+H_{1}\psi_{1}+H_{2}\psi_{0}=E_{0}\psi_{2}+E_{1}\psi_{1}+E_{2}\psi_{0}=E_{0}\psi_{2}+E_{2}\psi_{0}$.
We take again the inner product with $\psi_{0}$ and use the self-adjointness
of $H_{0}$ to obtain:

$\displaystyle
E_{2}=\frac{\langle\psi_{0},H_{1}\psi_{1}+H_{2}\psi_{0}\rangle}{\langle\psi_{0},\psi_{0}\rangle}=\langle\psi_{0},H_{1}\psi_{1}+H_{2}\psi_{0}\rangle.$

We now calculate $\psi_{1}$. It satisfies the equation

$\displaystyle H_{0}\psi_{1}+H_{1}\psi_{0}=E_{0}\psi_{1},$

which takes the form

$\displaystyle\psi_{1}^{\prime\prime}+m^{2}\psi_{1}=\frac{1}{\sqrt{\pi}}\cos(x)\cos(mx),$

in $(0,2\pi)$ with periodic boundary conditions. This is a second order
inhomogeneous ODE; its solution is:

$\displaystyle\psi_{1}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{(4m^{2}-1)\sqrt{\pi}}\Big{(}2m\sin(x)\sin(mx)+\cos(x)\cos(mx)\Big{)}+C_{2}\sin(mx)+C_{1}\cos(mx),$

for $C_{1},C_{2}\in\mathbb{R}$.

We can now calculate $\langle\psi_{0},H_{1}\psi_{1}\rangle$:

$\displaystyle\langle\psi_{0},H_{1}\psi_{1}\rangle=\begin{cases}-\frac{5}{12}\text{
if }m=1\\\ -\frac{1}{2(4m^{2}-1)}\text{ if }m\neq 1.\end{cases}$

Similarly, we have:

$\displaystyle\langle\psi_{0},H_{2}\psi_{0}\rangle=\begin{cases}-\frac{1}{2}\text{
if }m\neq 1\\\ -\frac{1}{4}\text{ if }m=1.\end{cases}$

We conclude that for $m=1$ the second order perturbation is
$E_{2}=-\frac{2}{3}$, and for $m\neq 1$:

$\displaystyle E_{2}=-\frac{1}{2(4m^{2}-1)}-\frac{1}{2},$

The eigenvalue of the perturbed operator $H$ is hence given by:
$E=E_{0}+\delta E_{1}+\delta^{2}E_{2}=-m^{2}+\delta^{2}E_{2}$.

We are primarily interested in the first nonzero eigenvalue. We plot the
asymptotic formula $E=-1-\frac{2}{3}\delta^{2}$ versus the numerically
obtained ones. The latter are obtained by using Matlab’s eig function on the
operator $Hf=f^{\prime\prime}-\delta\cos(x)f-\delta^{2}\sin^{2}(x)f$

Figure 2: First eigenvalue as a function of $\delta$, $W(x)=-\cos(x)$.

###  4.1 A higher order harmonic potential

We now consider the interaction potential $W(x)=-\cos(nx)$. The Schrödinger
operator becomes:

$\displaystyle
Hf=f^{\prime\prime}-\frac{n^{2}}{2}\left(\frac{\beta^{2}}{2}r^{2}\sin^{2}(nx)+\beta
r\cos(nx)\right)f(x).$

Writing $\beta=2(1+\varepsilon)$ and using the same expansion as before, we
obtain:

$\displaystyle Hf=f^{\prime\prime}-\delta
n^{2}\cos(nx)f(x)-\delta^{2}n^{2}\sin^{2}(nx)f(x)+O(\delta^{3})=:H_{0}+\delta
H_{1}+\delta^{2}H_{2}+O(\delta^{3}).$

We note that we can indeed use the same expansion for $r$ as before, as this
comes from using the Taylor expansion on $\exp(\beta r\cos(kx))$ and all terms
involving the $\cos(kx)$ terms simplify. More precisely, we now have:

$\displaystyle r=\frac{\int_{0}^{2\pi}\cos(nx)\exp(\beta
r\cos(nx))\,\mathrm{d}x}{\int_{0}^{2\pi}\exp(\beta r\cos(nx))\,\mathrm{d}x}.$

Approximating $e^{x}\approx 1+x+\frac{x^{2}}{2}$:

$\displaystyle r$ $\displaystyle=\frac{\int_{0}^{2\pi}\cos(nx)+\beta
r\cos(nx)^{2}+\frac{\beta^{2}}{2}r^{2}\cos(nx)^{3}\,\mathrm{d}x}{\int_{0}^{2\pi}1+\beta
r\cos(nx)+\frac{\beta^{2}}{2}r^{2}\cos(nx)^{2}\,\mathrm{d}x}=\frac{\beta
r\pi}{2\pi+\frac{\beta^{2}}{2}r^{2}\pi}.$

Substituting $\beta=2(1+\varepsilon)$, we deduce that

$\displaystyle r=\frac{\sqrt{\varepsilon}}{1+\varepsilon}.$

We study again the equations corresponding to different orders of
$\delta=\sqrt{\varepsilon}$.

Order $O(1)$: $H_{0}\psi_{0}=E_{0}\psi_{0}$. As before, this is just the
eigenvalue equation for the unperturbed operator $H_{0}f=f^{\prime\prime}$.
Therefore, $\psi_{0}(x)=-\frac{1}{\sqrt{\pi}}\cos(mx)$ for $m\in\mathbb{N}$.

Order $O(\delta)$: $H_{0}\psi_{1}+H_{1}\psi_{0}=E_{0}\psi_{1}+E_{1}\psi_{0}$.
Using the same reasoning as with the previous case, we obtain:

$\displaystyle
E_{1}=\langle\psi_{0},H_{1}\psi_{0}\rangle=\begin{cases}0&\text{ if }n\neq
2m\\\ -\frac{n^{2}}{2}&\text{ if }n=2m.\end{cases}$

Therefore, when $n$ is even, we now have a non-zero first order perturbation;
the corresponding eigenvalues of $H$ are:

$\displaystyle E=-\frac{n^{2}}{2}\left(\frac{1}{2}+\delta\right).$

Order $O(\delta^{2})$:
$H_{0}\psi_{2}+H_{1}\psi_{1}+H_{2}\psi_{0}=E_{0}\psi_{2}+E_{1}\psi_{1}+E_{2}\psi_{0}$.

As before:

$\displaystyle E_{2}=\langle\psi_{0},H_{1}\psi_{1}+H_{2}\psi_{0}\rangle$

We only look at the second perturbation in the case where $E_{1}=0$, i.e. when
$n\neq 2m$. Performing the same calculations as in the previous subsection, we
obtain that for $n\neq m$:

$\displaystyle
E_{2}=\frac{n^{2}}{(n^{2}-4m^{2})}\left(\frac{n^{2}}{2}\right)-\frac{n^{2}}{2},$

while for $n=m$:

$\displaystyle E_{2}$
$\displaystyle=-\frac{5n^{2}}{12}-\frac{n^{2}}{4}=\frac{2n^{2}}{3}.$

Therefore, we obtain the following eigenvalues for $H$:

  * •

For $n\in\mathbb{N}$ even,
$E_{n}=-\frac{n^{2}}{2}\left(\frac{1}{2}+\delta\right)$;

  * •

For $n,m\in\mathbb{N}$, $m\neq\frac{n}{2}$,
$E_{n}=-m^{2}+\delta^{2}\left(\frac{n^{2}}{(n^{2}-4m^{2})}\left(\frac{n^{2}}{2}\right)-\frac{n^{2}}{2}\right)$;

  * •

For $n\in\mathbb{N}$, $E_{n}=-n^{2}-\delta^{2}\left(\frac{n^{2}}{6}\right)$.

Table 1: Eigenvalues for $W(x)=-\cos(2x)$, $\delta=0.1$. m | Perturbation eigenvalue | Numerical eigenvalue |
---|---|---|---
1 | -1.2 | -1.2154 |
2 | -4.0267 | -4.0267 |
3 | -9.0225 | -9.0221 |
4 | -16.0213 | -16.0213 |
Table 2: Eigenvalues for $W(x)=-\cos(3x)$, $\delta=0.1$. m | Perturbation eigenvalue | Numerical eigenvalue |
---|---|---|---
1 | -0.964 | -0.9653 |
2 | -4.1029 | -4.1015 |
3 | -9.06 | -9.0600 |
4 | -16.0524 | -16.0524 |
Table 3: Eigenvalues for $W(x)=-\cos(4x)$, $\delta=0.1$. m | Perturbation eigenvalue | Numerical eigenvalue |
---|---|---|---
1 | -0.9733 | -0.9739 |
2 | -4.8 | -4.8615 |
3 | -9.144 | -9.1434 |
4 | -16.067 | -16.067 |

###  4.2 Multichromatic potentials

We now consider $W(x)=-\cos(x)-\frac{1}{2}\cos(2x)$. We recall that now the
stationary distribution is given by:

$\rho_{\infty}(x)=\frac{\exp(\beta(r_{1}\cos(x)+\frac{r_{2}}{2}\cos(2x)))}{\int_{0}^{2\pi}\exp(\beta(r_{1}\cos(x)+\frac{r_{2}}{2}\cos(2x)))\,\mathrm{d}x}$
(16)

with $r_{1}$ and $r_{2}$ satisfying the self-consistency equations:

$r_{1}=\frac{\int_{0}^{2\pi}\cos(x)\exp(\beta(r_{1}\cos(x)+\frac{r_{2}}{2}\cos(2x)))\,\mathrm{d}x}{\int_{0}^{2\pi}\exp(\beta(r_{1}\cos(x)+\frac{r_{2}}{2}\cos(2x)))\,\mathrm{d}x}$
(17)

and

$r_{2}=\frac{\int_{0}^{2\pi}\cos(2x)\exp(\beta(r_{1}\cos(x)+\frac{r_{2}}{2}\cos(2x)))\,\mathrm{d}x}{\int_{0}^{2\pi}\exp(\beta(r_{1}\cos(x)+\frac{r_{2}}{2}\cos(2x)))\,\mathrm{d}x}$
(18)

We need to find approximations for $r_{1}$ and $r_{2}$ near the critical value
of $\beta$. Using that $\exp(x)\sim 1+x+\frac{x^{2}}{2}$ in equations
(17)-(18), we obtain that, close to $\beta=2$, $r_{2}^{2}(\beta)$ behaves
linearly, and that

$r_{1}=\sqrt{r_{2}\left(r_{2}+\frac{2}{\beta}\right)}.$ (19)

To obtain an explicit expression for $r_{2}(\beta)$ we use a basic fitting
algorithm to obtain that the best approximation is given by

$\displaystyle r_{2}(\beta)=\sqrt{\frac{3}{2}-\frac{3}{\beta}}.$

Substituting this in (19), we obtain the corresponding approximating function
for $r_{1}$:

$\displaystyle
r_{1}(\beta)=\sqrt{\frac{3}{2}-\frac{3}{\beta}+\frac{2}{\beta}\sqrt{\frac{3}{2}-\frac{3}{\beta}}}$

We now use these approximations to study the Schrödinger operator for
$W=-\cos(x)-\frac{1}{2}\cos(2x)$. With $\rho_{\infty}$ as in (16), we have
that:

$\displaystyle
U(x):=\frac{\beta}{2}W*\rho_{\infty}(x)=-\frac{\beta}{2}r_{1}\cos(x)-\frac{\beta}{4}r_{2}\cos(2x)$

Therefore:

$\displaystyle\mathcal{L}f=f^{\prime\prime}-\left(\frac{(U^{\prime})^{2}}{4}+\frac{U^{\prime\prime}}{2}\right)f$
$\displaystyle=f^{\prime\prime}-\left(\frac{\beta^{2}}{4}(r_{1}\sin(x)+r_{2}\sin(2x))^{2}+\frac{\beta}{2}(r_{1}\cos(x)+2r_{2}\cos(2x))\right)f$

We now substitute $\beta=2(1+\varepsilon)$, $\eta={\varepsilon}^{1/4}$, use
the above expressions for $r_{1}(\beta),r_{2}(\beta)$ and ignore terms of
order higher than $O(\eta^{3})$. After long calculations, we obtain the
following expression for the perturbed operator:

$\displaystyle Hf$
$\displaystyle:=f^{\prime\prime}-\left(\sqrt{\frac{3}{2}}\eta^{2}\sin^{2}(x)+\left(\frac{3}{2}\right)^{1/4}\eta\cos(x)+2\sqrt{\frac{3}{2}}\eta^{2}\cos(2x)\right)f$
$\displaystyle=:H_{0}+\eta H_{1}+\eta^{2}H_{2},$

where

$\displaystyle H_{0}f=f^{\prime\prime};$ $\displaystyle
H_{1}f=-\left(\frac{3}{2}\right)^{1/4}\cos(x)f(x);$ $\displaystyle
H_{2}f=-\sqrt{\frac{3}{2}}\left(\sin^{2}(x)+2\cos(2x)\right)f(x)=-\frac{1}{2}\sqrt{\frac{3}{2}}(1+3\cos(2x))f(x).$

We do the usual perturbation argument to find the eigenvalues $E$ and
eigenfunctions $\psi$ solving $H\psi=E\psi$. We solve the equation order by
order.

Order $O(1)$: $H_{0}\psi_{0}=E_{0}\psi_{0}$.

As usual, this gives us $\psi_{0}=\frac{1}{\sqrt{\pi}}\cos(nx)$,
$E_{0}=-n^{2}$.

Order $O(\eta)$: $H_{0}\psi_{1}+H_{1}\psi_{0}=E_{0}\psi_{1}+E_{1}\psi_{0}$.
Performing the same steps as before gives us $E_{1}=0$ again.

Order $O(\eta^{2})$:
$H_{0}\psi_{2}+H_{1}\psi_{1}+H_{2}\psi_{0}=E_{0}\psi_{2}+E_{1}\psi_{1}+E_{2}\psi_{0}=E_{0}\psi_{2}+E_{2}\psi_{0}$.
We again take the inner product with $\psi_{0}$ and use that $H_{0}$ is
Hermitian to obtain:

$\displaystyle E_{2}=\langle\psi_{0},H_{1}\psi_{1}+H_{2}\psi_{0}\rangle$

Solving similar equations as in the previous subsections we have, for $n=1$:

$\displaystyle E_{2}=-\frac{(36)^{1/4}}{12}-\frac{3}{2}\sqrt{\frac{3}{2}}$

and for $n\neq 1$:

$\displaystyle
E_{2}=-\frac{1}{2}\sqrt{\frac{3}{2}}\left(\frac{1}{4n^{2}-1}+1\right)$

Therefore, the eigenvalues of $H$ are
$E=-1-\eta^{2}\left(\frac{(36)^{1/4}}{12}+\frac{3}{2}\sqrt{\frac{3}{2}}\right)$
and, for $n\neq 1$:

$\displaystyle
E=-n^{2}-\eta^{2}\frac{1}{2}\sqrt{\frac{3}{2}}\left(\frac{1}{4n^{2}-1}+1\right)$

Fixing $\eta$ (for example, $\eta=0.5$) and calculating a few values, we can
see that they are in very good agreement with the eigenvalues found
numerically by MATLAB. As before, we provide a plot to illustrate the
behaviour of the first eigenvalue as $\eta$ varies.

Figure 3: First eigenvalue as a function of $\eta$,
$W(x)=-\cos(x)-\frac{1}{2}\cos(2x)$.

##  5 Numerical experiments

In this section we consider both steady states and the dynamics of the Fokker-
Planck equation (1) by solving it numerically using pseudospectral methods [7,
28, 21]. Furthermore, we compare the solution of the PDE to results obtained
from Monte-Carlo simulations using the Euler-Maruyama method [23].

###  5.1 Steady states and intermediate dynamics

As described in Section 2, we expect there to be a qualitative change in the
nature of the steady-state solution of the PDE (1) at $\beta_{c}=T_{c}^{-1}$.
Due to the gradient structure of the PDE, (numerical approximations to the)
steady states can be determined either by solving the PDE directly over a long
time interval, or by solving the self-consistency equation (4) iteratively,
for example via Picard iteration. A key difference here is that iterative
approaches can converge to unstable steady states, whereas the PDE method
always approaches a stable steady state. This motivates our choice in this
section to use the long time PDE solutions.

In order to focus on the effects of the choice of the interaction parameters
on the dynamics, in this section we consider dynamics with no external
potential. We will reintroduce the external potential in the following section
when comparing against stochastic dynamics.

We reiterate from above that for for $\beta<\beta_{c}$, we expect the uniform
state to be stable, corresponding to the long time solution of the PDE being
$\rho_{\infty}=1/(2\pi)$. In contrast, for $\beta>\beta_{c}$, we expect to
observe the other, peaked solutions. Furthermore, as $\beta$ grows larger, we
expect the steady state to be more strongly peaked.

In Figures 4–7 we show ‘intermediate’ and ‘long’ time dynamics of the PDE for
various interactions, $W$, and a range of $\beta$. The particular choices of
$W$ are somewhat arbitrary; we have chosen interactions with $n=2$ and $n=3$
coefficients, with the values chosen to result in a range of different
dynamics.

We start with an initial condition that is a small perturbation of the uniform
state - $\rho_{0}(x)=\frac{1}{2\pi}+0.01\sin\left(x-\frac{\pi}{2}\right)$. We
note that the particular form of the initial condition dictates the position
of the peak in the steady state; we have chosen it such that the peak is
located at $\theta=\pi$ for clarity. Due to the translational invariance of
the problem, this perturbation can, in principle, be chosen to control the
position of the peak(s); we will make use of this in the following section
when comparing to the SDE dynamics. One may think that, in parameter regimes
where the uniform state is unstable, a deliberate perturbation is not
necessary to see dynamics, as this should be produced by the accumulation of
numerical errors. However, due to the very high accuracy of the chosen
numerical schemes we find that it is necessary to deliberately perturb the
uniform state in order to induce dynamics in reasonable computational times.

As can be seen in the figures, for $\beta<\beta_{c}$, the long time solution
is indeed the numerical steady state. As expected, for $\beta>\beta_{c}$ we
see a more interesting range of long time solutions, which (for the examples
below) have a single peak. Note that for $\beta$ close to $\beta_{c}$ the
convergence to equilibrium can be very slow; this is demonstrated in the right
hand plots of Figures 6 and 7 where the solutions for $\beta$ close to
$\beta_{c}$ have not yet converged to the final steady state. This effect is
due to the exponential slowing down of the dynamics near the critical value of
$\beta$.

Considering now the ‘intermediate’ dynamics in the left hand plots of Figures
6 and 7, we note that there are transient regimes where where the number of
peaks is the value of $k$ corresponding to the largest $|a_{k}|$. However,
these states appear to be unstable as the long time dynamics results in a
single peak. The instability of multi-peak steady states was also observed in
[33], [20].

Looking at the dynamics in a bit more detail, we see that for the interactions
in Figures 4 and 5, the largest (in magnitude) coefficient is $a_{1}$, which
corresponds to a single Fourier mode. As expected, this leads to relatively
simple dynamics in which the solution quickly converges to a single peak. As
shown in Figure 5 adding a third non-zero coefficient does not change the
critical value $\beta_{c}$, nor the qualitative nature of the solution.
However, it does affect the quantitative dynamics, for example by decreasing
the width of the final solutions for fixed $\beta$. There is also a more
obvious three-peaked state during the dynamics for $\beta=2$, although it is
possible that this arises at a different time for the dynamics in Figure 4.

In Figures 6 and 7 we consider interactions where the largest magnitude
coefficient corresponds to a higher Fourier mode. In these cases, we clearly
see the transient states which have a number of peaks corresponding to the
Fourier mode with the coefficient with the largest magnitude. However, the
steady state solution is dominated by the lowest non-zero Fourier mode, which
in the cases shown here results in a single peak.

Figure 4: PDE dynamics for different values of $\beta$ for the interacting
potential $W(\theta)=-3\cos(\theta)-\cos(2\theta),V(\theta)=0$ for times $T=5$
(left) and $T=10^{3}$ (right). The critical temperature is
$\beta_{c}=\frac{2}{3}$.

Figure 5: PDE dynamics for different values of $\beta$ for the interacting
potential $W(\theta)=-3\cos(\theta)-\cos(2\theta)-\cos(3\theta),V(\theta)=0$
for times $T=5$ (left) and $T=10^{3}$ (right). The critical temperature is
$\beta_{c}=\frac{2}{3}$.

Figure 6: PDE dynamics for different values of $\beta$ for the interacting
potential $W(x)=-\cos(x)-2\cos(2x),V(x)=0$, for times $T=10^{3}$ (left) and
$T=10^{6}$ (right). The critical temperature is $\beta_{c}=1$.

Figure 7: PDE dynamics for different values of $\beta$ for the interacting
potential $W(x)=-\cos(x)-\cos(2x)-3\cos(3x),V(x)=0$, for times $T=10^{3}$
(left) and $T=10^{6}$ (right). The critical temperature is
$\beta_{c}=\frac{2}{3}$.

Figures 4-7 only show stable steady states with one single peak. It is
possible to construct models with stable multipeak stationary distributions by
simply rescaling the domain of the interaction potential, i.e. by considering
potentials with higher harmonics and with a zero first Fourier mode. For
example, the interaction potential $W(x)=-\cos(2x)$ will have a two-peak
stable steady state.

###  5.2 Comparison with Monte Carlo Simulations

In this section we compare the dynamics of the PDE, (1), with those of the
underlying interacting particle SDE, (5). The aim here is twofold: we will (i)
demonstrate the effects of translational invariance on the results of
stochastic sampling, and (ii) compare the PDE and SDE dynamics directly for a
single initial condition with a large number of particles; this is a numerical
demonstration of the mean-field limiting dynamics agreeing.

We will also investigate systems which have a non-zero confining potential.
For a non-trivial potential, this breaks the translational invariance and
allows us to demonstrate very good agreement between the two solutions by
performing multiple runs of the stochastic dynamics and averaging. This is
computationally much cheaper than performing a single run with a larger number
of particles: For $N$ particles and $R$ runs, the computational cost scales
approximately as $N^{2}R$ (where the $N^{2}$ scaling results from having to
compute the interaction potential), so increasing the number of runs is much
cheaper than increasing the number of particles.

To numerically solve the SDE (5), we use the Euler-Maruyama method [23] and
compare our findings with the results obtained from the PDE solver. Unless
otherwise stated, in all the simulations we use 500 particles, and a timestep
of $0.01$. The initial condition is generated by sampling from the same
initial condition $\rho_{0}$ as for the PDE given in Section 5.1 using Monte-
Carlo slice sampling [27].

We first consider the interaction potential
$W(x)=-\cos(x)-\frac{1}{2}\cos(2x)$, with no confining potential, an
interaction strength $K=1$, and inverse temperature $\beta=3$. We run the
dynamics up to a final time $T=200$.

As mentioned, in the absence of a confining potential, the problem is
translationally invariant. This means, in particular, that above the phase
transition we have infinitely many stationary states parameterized by an angle
$\theta\in[0,2\pi]$ [5]. Consequently, when performing particle simulations,
if we average over many realizations of the noise, we obtain (approximately)
the uniform distribution; see the left plot in Figure 8. On the other hand, if
we perform a single simulation, with a sufficiently large number of particles,
then the results of the stochastic simulations are in good agreement with the
results of the PDE, up to a translation in space; see the right plot in Figure
8. Note that, in order to more clearly demonstrate the agreement, we have
adjusted the PDE initial condition through a translational shift so that the
positions of the peaks (approximately) align. This gives exactly the same
result as simply re-plotting the PDE (or SDE) solution on a translated axis.

Figure 8: PDE and SDE dynamics for $W(x)=-\cos(x)-\frac{1}{2}\cos(2x)$ for 500
particles and 100 runs (left) and 1 run (right).

We now compare the PDE and SDE dynamics for a range of other interactions and
no external potential. Our next example concerns the interaction potential
$W(x)=-\frac{1}{4}\cos(4x)-\frac{1}{6}\cos(6x)$, and is chosen to demonstrate
the agreement for multi-peaked solutions. We integrate up to $T=200$, with a
timestep $\Delta t=0.01$. As can be seen in Figure 9, the agreement between
the PDE and SDE is very good. Note that here we did not need to shift the PDE
solution to align the peaks; this would not be true for a different
realisation of the noise.

Figure 9: PDE and SDE dynamics with one run for
$W(x)=-\frac{1}{4}\cos(4x)-\frac{1}{6}\cos(6x),V(x)=0,\beta=10$.

Our next example is chosen from [35] concerning Hodgkin-Huxley oscillators,
where the second Fourier mode now has positive sign. We take
$W(x)=-\cos(x)+\frac{1}{2}\cos(2x)$. Our numerical experiments (not shown)
agree with the result of the paper that there is a phase transition at
$\beta_{c}=2$. We plot here the density for $\beta=3$, with $100$, T = $200$,
and $\Delta t=0.01$. The results are shown in Figure 10; again, the agreement
between the SDE and PDE is very good.

Figure 10: PDE and SDE dynamics for
$W(x)=-\cos(x)+\frac{1}{2}\cos(2x),V(x)=0,\beta=3$.

In contrast to the examples presented so far, the presence of an
external/confining potential breaks translation invariance. Consequently,
averaging over many realization of the noise in the particle simulations leads
to a non-uniform stationary state. We now consider two such examples,
averaging solutions for 500 particles over 10 runs in each case.

For the first example we consider the Brownian mean field model in a magnetic
field, (7), with $W(x)=-\cos(x)$ and $V(x)=0.2\cos(x+\pi)$. We integrate up to
$T=1000$ with $\Delta t=0.01$. The results are shown in Figure 11, where we
once again have excellent agreement between the PDE and SDE.

Figure 11: PDE and SDE dynamics for $W(x)=-\cos(x)$ with confining potential
$V(x)=0.2\cos(x+\pi),\beta=3$, 500 particles, 10 runs.

Finally for long time solutions, we consider an example taken from [1]. Here
$W(x)=-2\cos(2x)$ and $V(x)=\cos(x)$. In this case the dynamics converge more
slowly and we take a final time of $T=7000$ and $\Delta t=0.01$. Figure 12
demonstrates the excellent agreement between the dynamics at long times.

Figure 12: PDE and SDE dynamics for $W(x)=-2\cos(x)$ with confining potenital
$V(x)=\cos(x+\pi),\beta=3,500$ particles.

Of course, one is often interested not only in long time dynamics, but also in
the full evolution of the density. In Figure 13 we compare the PDE and SDE
dynamics at a number of times for two different combinations of interaction
and external potentials. As previously discussed, the agreement between the
two dynamics is better in the second picture as we have a confining potential
which breaks translation invariance, and a simpler interaction potential.

Figure 13: PDE and SDE dynamics for a different range of times, $\beta=3$, for
$W(x)=-\cos(x)-\frac{1}{2}\cos(2x),V(x)=0$ (left) and for $W(x)=-\cos(x)$,
$V(x)=0.2\cos\left(x-\frac{\pi}{2}\right)$ (right).

###  5.3 PDE simulations for the HKB model

Another model with confining potential is the Haken–Kelso–Bunz (HKB) (see
Chapter 5, [18]). Here $W(x)=-\cos(x)$,
$V(x)=-\alpha\cos(x+\pi)-\gamma\cos(2x)$. We plot the stationary states of the
Fokker-Planck PDE for different values of $\alpha,\gamma$.

Figure 14: PDE dynamics for the HKB model with $\alpha,\gamma=1$, $\beta=3$.
Figure 15: PDE dynamics for the HKB model with $\alpha=1,\gamma=\frac{1}{8}$,
$\beta=3$.

##  6 Conclusions

In this paper, we studied the stability of steady states for McKean-Vlasov
dynamics on the torus, for interaction potentials that contain multiple non-
zero Fourier modes. It was shown, through the study of the self-consistency
equations, of the free energy of the system and through a linear stability
analysis, that there is a critical temperature value at which the constant
stationary state of the system becomes unstable. Furthermore, the stability of
non-constant stationary states was analysed by means of perturbation theory
for higher harmonic and multichromatic interaction potentials. Finally, we
verified these results with extensive numerical simulations of both the
Fokker-Planck PDEs and the systems of SDEs involved.

The work presented in this paper can be extended in several interesting
directions. First, we would like to study the effect of inertia on the
formation and stability of multipeak solutions by considering the kinetic mean
field PDE. Second, the impact of colored noise on the stability on the phase
transitions and on the stability of non-uniform steady states is an important
question, motivated by recent work on the modeling of collective organization
in cyanobacteria [15]. More generally, we aim at applying our analytical and
numerical methodologies to the study of active matter, e.g. [32]. All these
topics are currently under investigation.

#### Acknowledgments

BB is funded by a studentship from the Imperial College London EPSRC DTP in
Mathematical Sciences (Grant No. EP/W523872/1) and was supported by a Mary
Lister McCammon summer research fellowship. GP is partially supported by an
ERC-EPSRC Frontier Research Guarantee through Grant No. EP/X038645, ERC
through Advanced Grant No. 247031.

## References

  * [1] L. Angeli, J. Barré, M. Kolodziejczyk and M. Ottobre  “Well-posedness and stationary solutions of McKean-Vlasov (S) PDEs”  In _Journal of Mathematical Analysis and Applications_ 526.2 Elsevier, 2023, pp. 127301
  * [2] R. Balescu  “Statistical dynamics. Matter out of equilibrium”  London: Imperial College Press, 1997
  * [3] G. A. Battle  “Phase transitions for a continuous system of classical particles in a box”  In _Comm. Math. Phys._ 55.3, 1977, pp. 299–315  DOI: 10.1007/BF01614553
  * [4] F. Bavaud  “Equilibrium properties of the Vlasov functional: the generalized Poisson-Boltzmann-Emden equation”  In _Rev. Modern Phys._ 63.1, 1991, pp. 129–148  DOI: 10.1103/RevModPhys.63.129
  * [5] L. Bertini, G. Giacomin and K. Pakdaman  “Dynamical aspects of mean field plane rotators and the Kuramoto model”  In _Journal of Statistical Physics_ 138 Springer, 2010, pp. 270–290
  * [6] J. Binney and S. Tremaine  “Galactic Dynamics”  Princeton: Princeton University Press, 2008
  * [7] J.P. Boyd  “Chebyshev and Fourier spectral methods”  Courier Corporation, 2001
  * [8] J. A. Carrillo, M. G. Delgadino and G. A. Pavliotis  “A $\lambda$-convexity based proof for the propagation of chaos for weakly interacting stochastic particles”  In _J. Funct. Anal._ 279.10, 2020, pp. 108734  DOI: 10.1016/j.jfa.2020.108734
  * [9] J. A. Carrillo, R. S. Gvalani, G. A. Pavliotis and A. Schlichting  “Long-time behaviour and phase transitions for the Mckean-Vlasov equation on the torus”  In _Arch. Ration. Mech. Anal._ 235.1, 2020, pp. 635–690  DOI: 10.1007/s00205-019-01430-4
  * [10] P. Chavanis  “The Brownian mean field model”  In _The European Physical Journal B_ 87 Springer, 2014, pp. 1–33
  * [11] P. Constantin, I. G. Kevrekidis and E. S. Titi  “Asymptotic states of a Smoluchowski equation”  In _Arch. Ration. Mech. Anal._ 174.3, 2004, pp. 365–384  DOI: 10.1007/s00205-004-0331-8
  * [12] F. A. Davidson and N. Dodds  “Spectral properties of non-local differential operators”  In _Appl. Anal._ 85.6-7, 2006, pp. 717–734  DOI: 10.1080/00036810600555171
  * [13] F. A. Davidson and N. Dodds  “Spectral properties of non-local uniformly-elliptic operators”  In _Electron. J. Differential Equations_ , 2006, pp. No. 126, 15
  * [14] M. G. Delgadino, R. S. Gvalani and G. A. Pavliotis  “On the Diffusive-Mean Field Limit for Weakly Interacting Diffusions Exhibiting Phase Transitions”  In _Arch. Ration. Mech. Anal._ 241.1, 2021, pp. 91–148  DOI: 10.1007/s00205-021-01648-1
  * [15] M. K. Faluweki, J. Cammann, M. G. Mazza and L. Goehring  “Active Spaghetti: Collective Organization in Cyanobacteria”  In _Phys. Rev. Lett._ 131 American Physical Society, 2023, pp. 158303  DOI: 10.1103/PhysRevLett.131.158303
  * [16] I. Fatkullin and V. Slastikov  “A note on the Onsager model of nematic phase transitions”  In _Commun. Math. Sci._ 3.1, 2005, pp. 21–26  URL: http://projecteuclid.org/euclid.cms/1111095638
  * [17] I. Fatkullin and V. Slastikov  “Critical points of the Onsager functional on a sphere”  In _Nonlinearity_ 18.6, 2005, pp. 2565–2580  DOI: 10.1088/0951-7715/18/6/008
  * [18] T.D. Frank  “Nonlinear Fokker-Planck equations”, Springer Series in Synergetics  Berlin: Springer-Verlag, 2005, pp. xii+404
  * [19] J. Garnier, G. Papanicolaou and T. Yang  “Consensus convergence with stochastic effects”  In _Vietnam Journal of Mathematics_ 45 Springer, 2017, pp. 51–75
  * [20] E. Geigant and M. Stoll  “Stability of peak solutions of a non-linear transport equation on the circle”  In _Electronic Journal of Differential Equations_ 2012.157 Citeseer, 2012, pp. 1–41
  * [21] B. D. Goddard, A. Nold and S. Kalliadasis  “2DChebClass [Software]”  Edinburgh DataShare, 2017  URL: http://dx.doi.org/10.7488/ds/1991
  * [22] B. D. Goddard, B. Gooding, H. Short and G. A. Pavliotis  “Noisy bounded confidence models for opinion dynamics: the effect of boundary conditions on phase transitions” hxab044  In _IMA Journal of Applied Mathematics_ , 2021  DOI: 10.1093/imamat/hxab044
  * [23] P.E. Kloeden and E. Platen  “Numerical Solution of Stochastic Differential Equations”, Stochastic Modelling and Applied Probability  Springer Berlin Heidelberg, 2013  URL: https://books.google.co.uk/books?id=r9r6CAAAQBAJ
  * [24] D. Lacker  “Hierarchies, entropy, and quantitative propagation of chaos for mean field diffusions”  In _Probab. Math. Phys._ 4.2, 2023, pp. 377–432  DOI: 10.2140/pmp.2023.4.377
  * [25] M. Lucia and J. Vukadinovic  “Exact multiplicity of nematic states for an Onsager model”  In _Nonlinearity_ 23.12, 2010, pp. 3157–3185  DOI: 10.1088/0951-7715/23/12/009
  * [26] N. Martzel and C. Aslangul  “Mean-field treatment of the many-body Fokker-Planck equation”  In _J. Phys. A_ 34.50, 2001, pp. 11225–11240  DOI: 10.1088/0305-4470/34/50/305
  * [27] R. M. Neal  “Slice sampling”  In _The annals of statistics_ 31.3 Institute of Mathematical Statistics, 2003, pp. 705–767
  * [28] A. Nold et al.  “Pseudospectral methods for density functional theory in bounded and unbounded domains”  In _J. Comput. Phys._ 334 Elsevier, 2017, pp. 639–664
  * [29] K. Oelschlager  “A martingale approach to the law of large numbers for weakly interacting stochastic processes”  In _The Annals of Probability_ JSTOR, 1984, pp. 458–479
  * [30] K. J. Painter, T. Hillen and J. R. Potts  “Biological modeling with nonlocal advection-diffusion equations”  In _Math. Models Methods Appl. Sci._ 34.1, 2024, pp. 57–107  DOI: 10.1142/S0218202524400025
  * [31] G. A. Pavliotis  “Stochastic processes and applications” Diffusion processes, the Fokker-Planck and Langevin equations 60, Texts in Applied Mathematics  Springer, New York, 2014, pp. xiv+339  DOI: 10.1007/978-1-4939-1323-7
  * [32] F. Peruani, A. Deutsch and M. Bär  “A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms”  In _The European Physical Journal Special Topics_ 157 Springer, 2008, pp. 111–122
  * [33] I. Primi, A. Stevens and J.L. Velazquez  “Mass-selection in alignment models with non-deterministic effects”  In _Communications in Partial Differential Equations_ 34.5 Taylor & Francis, 2009, pp. 419–456
  * [34] J. Vukadinovic  “Inertial manifolds for a Smoluchowski equation on the unit sphere”  In _Comm. Math. Phys._ 285.3, 2009, pp. 975–990  DOI: 10.1007/s00220-008-0460-2
  * [35] J. Vukadinovic  “Phase transition for the McKean-Vlasov equation of weakly coupled Hodgkin-Huxley oscillators”  In _Discrete and Continuous Dynamical Systems_ DiscreteContinuous Dynamical Systems, 2023, pp. 0–0


real spacetime manifold $\mathcal{M}$ provides a dynamical dimension not
necessarily equal to the standard (Lebesgue covering) dimension
$n=\dim\mathcal{M}.$ So consider the $n$-dimensional diffusion equation
$\partial_{s}\,\phi(s,x)=-D\,\Box^{(n)}(x)\,\phi(s,x),$ where $s$ is some
fictitious (nonphysical) “time,” $\Box^{(n)}\,\phi\doteq P[\phi]$ is the
$n$-dimensional massless scalar wave operator on $\mathcal{M}$ such as eqn.
(184), and $D$ is a positive real (diffusion) constant. The initial condition
at $s=0$ is $\phi(0,x)=\delta^{(n)}(x).$ On standard $n$-dimensional Euclidean
space, $\Box^{(n)}$ is the $n$-dimensional Laplacian, and one has the well-
known the Gaussian solution

$\displaystyle\phi(s,x)$
$\displaystyle=\Big{(}\frac{1}{2\pi\sigma(s)^{2}}\Big{)}^{n/2}\exp\big{[}-x^{2}/2\sigma(s)^{2}\big{]},\quad\text{and}$
(188) $\displaystyle\sigma(s)^{2}$ $\displaystyle=2Ds.$ (189)

If $\phi(s,x)$ is interpreted as the probability of return to $x$ by “time”
$s,$ then the spectral dimension $n_{s}$ can be read off from (188) as

$n_{s}=-2\lim_{s\to+\infty}\ln\phi(s,0)/\ln s\,.$ (190)

Since the initial simulations demonstrated spectral dimension reduction,
researchers have provided phenomenological insights into this intriguing
behavior [32]. In particular, interpretations involving the large momentum
scaling of the scalar wave operator $P[\phi]$ have been presented. So far all
these descriptions are quantum in nature. For classical DGR there is a
deformation of the usual scalar wave operator due to its internal
$\star$-products. According to the IKT representation of $\star$ there are
special $N_{X}$-dimensional submanifolds, namely the TSM. Could the the
averaging _only_ on the TSM lead to classical spectral dimension reduction on
the length scale $\sim\xi_{c}\gg L_{P}?$

The massless $\star$-scalar wave operator of interest here was derived earlier
as (184), namely

$P[\phi]=\big{[}g^{\mu\nu}\star\partial_{\nu}\,\phi\big{]}_{\stackrel{{\scriptstyle\star}}{{,}}\mu}$
(191)

in any coordinate system. One is interested in the scaling of this operator at
large $4$-momentum $k^{\mu}$ of $\phi.$ Since the $\star$-inverse metric
$g^{\mu\nu}$ is roughly like the reciprocal of an area $A,$ one can crudely
conceive of this as a scaling of $1/A(\ell)$ with length $\ell.$ In the case
of quantum gravity, one makes use of the loop quantum approach result that the
area operator has minimum eigenvalues of order $L_{P}^{2},$ and sets
$A(\ell)=\ell^{2}+L_{P}^{2}.$ Putting $k\simeq 1/\ell,$ the inverse-area
scaling is described by

$\mathcal{F}(k)=\big{[}A(1/k,L_{P})\big{]}^{-1}/\big{[}A(1/k_{0},L_{P})\big{]}^{-1}$
(192)

where $k_{0}\ll L_{P}$ is some “small” fiducial momentum where one fixes
$\mathcal{F}(k_{0})=1.$ Then

$\mathcal{F}(k)=\frac{(k^{-2}+L_{P}^{2})^{-1}}{(k_{0}^{-2}+L_{P}^{2})^{-1}}=\frac{k^{2}}{1+(L_{P}\,k)^{2}}\cdot\frac{1+(L_{P}\,k_{0})^{2}}{k_{0}^{2}}.$
(193)

One finds $\mathcal{F}(k)\simeq(k/k_{0})^{2}$ for $k_{0}\ll k\ll 1/L_{P}$ and
$\mathcal{F}\simeq(k_{0}L_{P})^{-2}$ for $k\gg 1/L_{P}\gg k_{0}.$ The
$k$-independent saturation of $\mathcal{F}(k)$ for $k\gtrsim L_{P}$ is the
hallmark of spectral dimension reduction at the Planck scale.

Here we apply this phenomenological approach to DGR. $\xi_{c}$ takes on the
role of $L_{P}$ because as the initial condition for $\star$-diffusion
$\delta^{\star}_{x}$ acts approximately as a Gaussian of extent $\xi_{c}$ only
on the TSM$(x).$ Also approximate
$P_{\phi}(k)\simeq-g^{\mu\nu}\star(k_{\mu}k_{\nu}\phi).$ Separating
$P_{\phi}(k)$ into its action on the $N_{X}$-dimensional TSM and the
orthogonal $(n-N_{X})$ directions, one may write

$\displaystyle\mathcal{F}(k)$
$\displaystyle\simeq\frac{1}{n}\,\Big{[}\;\sum_{a=1}^{N_{X}}\frac{k_{a}^{2}}{1+(\xi_{c}\,k_{\parallel})^{2}}\frac{1+(\xi_{c}\,k_{0})^{2}}{k_{0}^{2}}+\sum_{a=N_{X}+1}^{n}\frac{k_{a}^{2}}{k_{0}^{2}}\;\Big{]}$
(194)
$\displaystyle\simeq\frac{1}{n}\,\Big{[}\,\frac{N_{X}\,k_{\parallel}^{2}}{1+(\xi_{c}\,k_{\parallel})^{2}}\frac{\big{(}1+(\xi_{c}\,k_{0})^{2}\big{)}}{k_{0}^{2}}+\frac{(n-N_{X})\,k_{\perp}^{2}}{k_{0}^{2}}\,\Big{]}.$
(195)

Here $k_{\parallel}$ refers to the projection of $k$ onto the TSM, $k_{\perp}$
is the off-TSM part. The overall $1/n$ factor maintains $\mathcal{F}(k_{0})=1$
for $k_{0}\ll 1/\xi_{c}.$ For $k_{0}\ll|k_{a}|\ll\xi_{c}^{-1},$
$\mathcal{F}(k)\simeq k^{2}/k_{0}^{2}.$ For $\phi$ propagation on the TSM one
drops the second term on the RHS of (195) and

$\mathcal{F}_{\parallel}(k)\simeq\frac{1}{n}\,\Big{[}\,\frac{N_{X}\,k_{\parallel}^{2}}{1+(\xi_{c}\,k_{\parallel})^{2}}\,\frac{1}{k_{0}^{2}}\,\Big{]}\to\frac{N_{X}}{n(\xi_{c}\,k_{0})^{2}}\quad\text{for}\quad
k_{\parallel}\gg\xi_{c}^{-1}.$ (196)

This large $k$ behavior saturates, signaling dimensional reduction. Similarly
off TSM one drops the first term to obtain

$\mathcal{F}_{\perp}(k)\simeq\frac{n-N_{X}}{n}\,\Big{(}\frac{k_{\perp}}{k_{0}}\Big{)}^{2},$
(197)

which never saturates, and so is not dimensionally reduced. For $k\simeq
k_{\perp}\simeq k_{\parallel}$ one gets

$\mathcal{F}(k\simeq k_{\perp}\simeq
k_{\parallel})\simeq\frac{n-N_{X}}{n}\,\Big{(}\frac{k}{k_{0}}\Big{)}^{2}\quad\text{for}\quad
k\gg\xi_{c}^{-1}$ (198)

a non-saturating limit and thus not dimensionally reduced. Hence one has found
classical anisotropic scalar diffusion, but not an overall reduction of
spectral dimension as in the quantum case. Only diffusion on the TSM exhibits
spectral dimension reduction. This micro-anisotropy of spacetime will be
explored further in the next section.

##  9\. Experimental Tests of DGR

Due to the non-locality of DGR both the celebrated Spin-Statistics Theorem
(SST) and CPT-symmetry from standard quantum field theory are expected to be
violated. As a result of the breakdown of the spin-statistics relationship,
one also anticipates deviations from the Pauli Exclusion Principle, for which
there are stringent experimental bounds [33]. One way to approach SST
breakdown is to study the parameter $\alpha_{ij}$ governing the exchange of
identical particles $i,j$ in a multi-particle state $\Phi.$ Here we assume
that spacetime is adiabatic in the Birrell-Davies sense [34]: The
characteristic frequencies $\omega$ of particles (matter modes) are much
larger than the rates of change of metrical variables, such as $H(a).$ It is
then permissible to speak of particles with their associated quantum numbers
$q$ that comprise the state $\Phi,$ which is taken to be a ray of vectors in
some Hilbert space.

Suppose one exchanges two identical particles $i$ and $j$ in $\Phi,$ so
$\Phi\to\alpha_{ij}\,\Phi$ for complex $\alpha_{ij}.$ If $\Phi$ is normalized
to unity, then $|\alpha_{i,j}|=1.$ Upon repeating the exchange, one returns to
the original state, hence $\alpha_{ij}\,\alpha_{ji}=1.$ In standard quantum
field theory $\alpha_{ij}$ is taken to be a function of a _single_ particle
property, namely to depend solely on the species $n=n_{i}=n_{j}$ of particles
$i$ and $j.$ This way $\alpha_{ij}(n)=\alpha_{ji}(n),$ and then
$(\alpha_{ij})^{2}=1,$ yielding $\alpha_{ij}=\pm 1$ and the familiar Bosonic
and Fermionic quantum statistics. In a non-local theory with characteristic
length $\xi_{c},$ $\alpha_{ij}$ may depend on other particles within a
distance $\sim\xi_{c}$ of $i,j$, as “influential spectators” $s.$ So
$\alpha_{ij}$ depends not only on the quantum numbers $q_{i},q_{j}$ of $i,j,$
respectively, but possibly also on $q_{s};$ i.e., on all the quantum numbers
specifying $\Phi.$ $\alpha$ is no longer simply a single particle property of
the exchanged particles. However if one arranges all those quantum numbers
into a list, then the ordering of that list of numbers on which $\alpha$
depends is not expected to be physically relevant. That is, although $\alpha$
itself has well-known physical effects, the value of $\alpha$ is not
anticipated to depend on how theorists choose to list the quantum numbers on
which $\alpha$ depends. Making that assumption, $\alpha_{ij}$ must depend on a
symmetrized combination of all the quantum numbers specifying $\Phi.$ In
particular, $\alpha_{ij}=\alpha_{ji},$ once again arriving at $\alpha_{ij}=\pm
1.$ Here one is also assuming that $\alpha$ does not depend on any paths taken
to assemble $\Phi,$ such as anyons which can only occur in two spatial
dimensions. As Weinberg [35] has pointed out, $\alpha_{ij}$ cannot depend only
on the spins $\sigma_{i},\sigma_{j},\sigma_{s}$ because then $\alpha$ would
form a $1$-dimensional representation of the rotation group, and there are
none. Nevertheless, there is still room for novel effects. Even with
symmetrized $\alpha$ in $3+1$-dimensions, one could have
$\alpha_{ij}=f((p_{i}-p_{j})^{2}),$ some function of the $4$-momentum
exchange. But since $\alpha_{ij}=\pm 1,$ any changes in $\alpha_{ij}$ must be
_discontinuous._ This discontinuous change would be expected to occur in a
non-local theory at characteristic $4$-momenta $\sim\xi_{c}^{-1}.$ The
possibility that $\alpha$ is no longer simply a single particle species
dependent constant causes multi-particle states to possess a potential “split-
personality disorder” as far as their statistics are concerned. A state
comprised of two identical spin-$0$ particles may behave quite decorously (Dr.
Jekyll) as 2 bosons below some threshold $4$-momentum exchange, only to act as
a nefarious fermion pair (Mr. Hyde) above threshold! Within DGR $\xi_{c}^{-1}$
is expected to lie in the range $(10^{-6}-10^{-2})E_{P},$ and experiments at
smaller characteristic energies will be censored from witnessing the
personality changes. Specifically, this would apply to the tests of the Pauli
Exclusion Principle performed to date, since their typical energies are on the
atomic to nuclear scales, at most $10-100$ MeV, many orders of magnitude
beneath any expected non-locality induced threshold(s).

The familiar ladder operator $\hat{a}(q),\hat{a}^{\dagger}(q)$ algebra for a
particle with quantum numbers $q$ is derived from the $\alpha_{ij}$ taken as a
single particle species property [36]. For instance one constructs
$N$-particle states from the particle vacuum state $\Phi_{0}$ as
$\hat{a}^{\dagger}(q_{1})\cdots\hat{a}^{\dagger}(q_{N})\Phi_{0}.$ In the
standard ladder algebra, where $\hat{a}(q)$ is defined as the adjoint of
$\hat{a}^{\dagger}(q),$ if $q,q_{1},\ldots,q_{N}$ are all sub-threshold bosons
(upper sign choice) or all fermions (lower sign) one finds

$\hat{a}(q)\;\Phi_{q_{1}\ldots q_{N}}=\sum_{r=1}^{N}\,(\pm
1)^{r+1}\,\delta(q-q_{r})\,\Phi_{q_{1}\ldots q_{r-1}q_{r+1}\ldots q_{N}}.$
(199)

Above threshold this becomes disrupted as soon as $\alpha_{ij}$ becomes
dependent on the influential spectators $s\neq i,j$ in an $N$-particle state,
since then $\alpha_{ij}$ can be $N$-dependent. Consequently above threshold
$\hat{a}(q)\,\Phi_{0}\neq 0$ becomes possible along with

$[\hat{a}(q^{\prime}),\hat{a}^{\dagger}(q)]_{\mp}\;\Phi_{q_{1}\ldots
q_{N}}\neq\delta(q^{\prime}-q)\,\Phi_{q_{1}\ldots q_{N}}$ (200)

and $[\hat{a}(q^{\prime}),\hat{a}(q)]_{\mp}\neq 0.$ The physical process
underlying this is that the ladder algebra encodes the cluster decomposition
principle. This is violated by DGR since
$[\hat{\psi}(x),\hat{\psi}(y)]_{\mp}\neq 0$ for space-like separated $x,y$
within a proper distance $\sim\xi_{c}$ of each each other. The ladder algebra
is an inappropriate foundation for quantum field theory above threshold. In
that strange land one must bid adieu to the familiar Feynman diagram
technology, loops, and so on.

While spin-statistics, the Pauli Exclusion Principle, and the ladder operator
algebra arise from particle exchange effects on multi-particle _states,_ CPT
symmetry is a symmetry of a local Lorentz invariant _action_ under
transformation of the _fields_ as $\psi\to(CPT)\,\psi.$ In a non-local theory
CPT symmetry will be broken. But unlike $\alpha_{ij}=\pm 1$ or
$|\alpha_{ij}|=1$ for particle exchange, there are no corresponding
thresholding restrictions for CPT violations. In fact it is easy to show that

$\displaystyle(CPT)\big{\langle}\hat{X}^{\mu}(x)\big{\rangle}(CPT)^{-1}$
$\displaystyle\neq\big{\langle}\hat{X}^{\mu}(x)\big{\rangle}\quad\text{and}$
(201) $\displaystyle(CPT)(\star_{X})(CPT)^{-1}$
$\displaystyle\neq(\star_{X}).$ (202)

What is expected size of the CPT violations? In an experiment with some
typical $4$-momentum $p,$ the first order $\star$-interactions relative to the
standard dot-product ones is

$\eta(\mathrm{CPT})\simeq\frac{\xi_{c}^{2}\,\|p\|}{L_{c}},$ (203)

with $L_{c}$ the gravitational curvature radius at the experiment. On the
Earth’s surface, $L_{c}\simeq(c^{2}R^{3}/GM)^{1/2}\simeq 10^{11}\,\mathrm{m},$
and then $\eta(CPT)\simeq(5\cdot
10^{-53})(\xi_{c}/10^{5}L_{P})^{2}(E(p)/\mathrm{TeV}).$ Near a black hole with
$L_{c}\simeq 10^{4}$ m, $\eta(CPT)$ is enhanced by a factor $10^{7},$ but even
that really does not help the experimental situation much. At the end of the
ACE in the very early Universe, $\eta(CPT)\sim 1.$

The current experimental bounds on $\eta(CPT)$ are too insensitive to detect
such minuscule violations. The Brookhaven measurement of the sidereal
variations of the muon magnetic moment’s $g-2$ found no detectable violations
down to $10^{-24}$ GeV, or $\eta(CPT)\lesssim 10^{-23}$ at $100$ MeV energies
[37]. The corresponding DGR estimate for flat spacetime is of second order in
momentum: $\eta\sim(\|p\|\,\xi_{c})^{2}\lesssim 10^{-28}.$ At larger energies,
experiments bounding the mass difference between the top and anti-top quark
masses to $67.1$ ppm at $7$ TeV imply $\eta(CPT)\lesssim 6.7\cdot 10^{-5}$ at
that energy [38]. Even if CPT violations are ever detected, they could arise
from violations of Lorentz invariance or non-locality, so by themselves they
would not constitute unambiguous evidence for DGR.

One may also consider deviations from maximum particle velocities being
bounded by $c$ that could arise from $\star$-wave equations and their
violation of the cluster decomposition principle. These can only occur in
curved spacetimes, coming from the kinetic term in the $\star$-action. This
results in the current experimental bounds being at least $24$ orders too
insensitive to the $\star$-corrections relative to the undeformed
contributions, which theoretical ratio is about the same as $\eta(CPT)$ from
eqn. (203).

Yet another class of possible experimental tests are modern Michelson-Morley
type experiments that search for anisotropies in the phase velocity of
photons. DGR can produce an effect like this but again only in curved
spacetimes as the photon propagation direction varies with respect to the
relative velocity of the local flow of GM particles. The 2009 experimental
bound [39] on this anisotropy is an impressive $\lesssim 10^{-17}.$ On the
other hand, the DGR estimate for photon wavelength $\lambda$ on Earth is

$\eta(MM)\simeq\frac{\xi_{c}^{2}}{L_{c}\lambda}\simeq(5.1\cdot
10^{-65})\,(\xi_{c}/10^{5}L_{P})^{2}\,(\lambda/500\,\mathrm{nm})^{-1},$ (204)

still a long way beyond present technology.

Hughes-Drever experiments are precision spectroscopic measurements of the
anisotropy of space and mass. First performed in 1960-1, they are among the
most accurate tests of relativity.[40] The twist might affect nuclear
transitions, here between Zeeman levels in an external magnetic field, because
the transition frequencies could potentially become anisotropic respect to the
twist generator(s). By comparing resonant Larmor frequencies in a
${}^{129}\mathrm{Xe}/{}^{3}\mathrm{He}$ co-magnetometer, Allmendinger et al.
[41] determined $\Delta M/M<6.7\cdot 10^{-34}$ for the nucleon. Both these
nuclei are spin $1/2,$ so by the Wigner-Eckart theorem they have no
electric/magnetic quadrupole or higher moments. A twist produced externally to
the nuclei would alter the nuclear magnetic moments so the ratio of their
Larmor frequencies in the same external magnetic field would vary with time as
$X$ in the lab frame changed with the Earth’s sidereal or annual motions. If
the twist arises from GM matter that is the Dark Matter comprising the
galactic halo, then $(\Delta
M/M)_{\mathrm{DGR}}\simeq(\xi_{c}(\mathrm{meas})/\Lambda_{C}(N))^{2}\simeq(\xi_{c}(\mathrm{meas})/10^{3}\,L_{P})^{2}(10^{-32}).$
Here $\Lambda_{C}(N)$ is the nucleon Compton wavelength, and
$\xi_{c}(\mathrm{meas})=\xi_{c}(X)\,(\beta\gamma),$ where the special
relativistic factor $(\beta\gamma)\simeq 6.9\cdot 10^{-4}$ enters because
$X=X^{0}\,\partial_{0}$ in the GM or halo’s rest frame while the Earth moves
through it along its galactic orbit at $v\simeq 206$ km/sec. Allmendinger et
al’s 2013 $\Delta M/M$ bound implies $\xi_{c}(X_{\mathrm{DM}})/L_{P}<3.76\cdot
10^{5},$ consistent with DGR expectations. On the other hand, if stable
Standard Model particles such as nucleons or nuclei generate the twist in
their own rest frame, then by symmetry that local twist generator vector must
lie parallel/anti-parallel to the particle’s spin, even during the spin’s
resonant dynamics. Hence in that case the magnetometer nuclei will perceive no
sidereal or annual shifts in the ratio of their Larmor precession frequencies,
producing a null experimental result, and no bound on $\xi_{c}(X).$
Nevertheless, future $4-5$ orders of magnitude improvements in Hughes-Drever
type experimental precision and accuracy could possibly rule out $X$ coming
from non-standard matter by bounding
$\xi_{c}(X_{\mathrm{DM}})/L_{P}\lesssim\mathscr{O}(1).$

##  10\. Conclusions and perspective

So where are we? What has been learned about DGR, what remains obscure? DGR
has been outfitted with a curved spacetime compatible integral kernel, whose
convergence and equivalence to an extended DDT have been clearly delineated.
Practical working details of $\star$-coord reparams ($\star 4$-diffs) have
been laid out and applied to develop firm foundations for classical
$\star$-gauge theories ($\star$-Cartan structure and $\star$-Bianchi identity)
and $\star$-Einstein field equations for simple dynamics $\delta
X=0=\delta(\star).$ The non-simple case was demonstrated to be non-viable as a
classical field theory. Photon and graviton ghosts were banished from
Minkowski spacetime. $\star$-FLRW cosmological models were studied, and an
explicit solution was found in the radiation dominated epoch for the spatially
flat case without a cosmological constant. BBN constraints were applied to
constrain the expressions for the twist generators $X.$ DGR does provide
viable models for dark matter as the twist producing GM matter such as the
Katyusha forms for $X.$ In those cases, dark matter would likely be a plasma
of dark-$U(1)$ interacting scalars or fermions. DGR can also qualitatively
account for the matter-antimatter asymmetry of the observable Universe. It
further predicts thresholded violations of the Spin-Statistics Theorem and
Pauli Exclusion Principle that lead to an unexpected Dr. Jekyll-Mr. Hyde high
energy statistical behavior. Future Hughes-Drever measurements hold out some
promise as experimental tests of whether twists are produced by standard
matter or not.

However for each of these “successes” there are corresponding unresolved
issues, no-shows, or loud calls of “out of bounds.” When applied to cosmology
DGR’s $\star$-FLRW models provide no resolution for the horizon or flatness
problems. There is no $\star$-ADM $3+1$ decomposition of $\star 4$-diffs, and
the use of gauge holonomies is restricted to loops. The elegant methods of
knot invariants and spin networks, so successful for quantizing canonical GR,
are unfortunately useless for DGR. Micro-causality and the cluster
decomposition principle are strongly violated on the $\xi_{c}(a=1)\simeq
10^{3-5}\,L_{P}$ length scale, yet attempts to explicitly calculate how those
cornerstones of quantum field theory are recovered at longer lengths were all
rebuffed. The intriguing non-canonical deformation flow process
$X^{(n)}_{\mathrm{sc}}\to X^{(n+1)}_{\mathrm{sc}}$ also remains a challenging
unsolved problem. Looking towards future theoretical developments, the method
of algebraic quantum field theory has shown some recent progress as being
adaptable to $\star$-manifolds with fixed classical $X_{\mathrm{sc}}.$ It is
able to enter into the (Birrell-Davies) non-adiabatic regime, where the
particle concept (ladder operator algebra) and micro-causality are not valid,
but regrettably it does not address deformation flow. A detailed understanding
of the predicted non-locality dominated early cosmological acausal epoch (ACE)
has posted “No Trespassing!” writ large for classical DGR. Finally, aside from
the potentially hopeful exceptions of Hughes-Drever and CPT-violation
measurements, the possible experimental tests of DGR examined so far lie many
orders of magnitude beyond present technology.

To place these virtues and defects of DGR into better perspective, we present
a brief discussion within the larger context of non-local modified theories of
gravity. These kinds of theories have been considered at least since the late
1980’s.[42][43] They consist of gravitational actions made from an infinite
set of higher derivative operators acting on curvature terms, and go by the
name of infinite derivative gravity (IDG) theories. Theories that just add a
finite number of higher order curvature terms like $R^{2}$ to the Einstein-
Hilbert action are necessarily local. The most general 4-dimensional, torsion-
free, and parity symmetric IDG action around a constant curvature manifold was
derived by [44], and reads

$S=\int\mathrm{d}^{4}x\,|g|^{1/2}\,\Big{[}M_{P}^{2}\,R+R\,F_{1}(\Box)\,R+R^{\mu\nu}\,F_{2}(\Box)\,R_{\mu\nu}+C^{\mu\nu\lambda\sigma}\,F_{3}(\Box)\,C_{\mu\nu\lambda\sigma}\,\Big{]},$
(205)

with $C$ the Weyl tensor, $F_{i}(\Box)=\sum_{k=0}^{\infty}f_{ik}(\Box/M^{2})$,
numerical coefficients $f_{ik}$, $\Box=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$
(d’Alembertian), and UV mass scale $M$. To be ghost free, the $F_{i}$ must all
be exponentials of entire functions.[45] Bounds on $M$ were obtained by [46]
and [47] as $M\gtrsim 10^{-5}\,E_{P}$, the same order of magnitude as
$\xi_{c}^{-1}$ in DGR. In contrast to IDG, theories have also been developed
possessing infra-red nonlocalities, partly to understand the cosmic
acceleration and/or cosmological constant problems.[48] There also have been
results exploring infra-red nonlocal covariant actions that include MOND
gravity.[49] Since DGR is UV nonlocal, we focus on comparing it to IDG. Both
theories are background independent and respect (local) Lorentz invariance.
However, IDG is invariant under standard local diffs, while DGR is invariant
under GM matter controlled, nonlocal $\star$-diffs. IDG is UV-complete and
(super) renormalizable/finite [50]; but classical DGR is an effective non-
renormalizable theory, valid only up to energy scales $\sim\xi_{c}^{-1}$.
IDG’s milder covariant d’Alembertian-based nonlocality leaves the light cone
invariant and preserves the cluster decomposition principle and micro-
causality. By contrast, DGR is more disruptive and violates these on length
scales $\lesssim\xi_{c}$. IDG’s nonlocality only affects the gravitational
sector, while DGR’s $\star$ enters all interactions in order for the total
action to be $\star$-diff invariant. The standard ADM $3+1$ decomposition is
expected to be valid for IDG since it has been reported to possess a
Hamiltonian formulation.[51] That decomposition does not obtain for DGR. So
far IDG has made no statements about dark matter, but DGR provides natural
models for it. In the large energy limit, IDG displays quantum spectral
dimension reduction, whereas DGR cannot reliably attain that quantum regime,
but instead shows classical anisotropic diffusion effects. Since DGR violates
the cluster decomposition principle on short length scale, it can accommodate
early universe matter-antimatter asymmetry, but IDG is silent about this. Both
IDG and DGR nonlocalities violate the Raychaudhuri null focussing theorem.
Nevertheless, IDG still exhibits spacetime solutions that include
singularities.[52] On the other hand, DGR encounters
predictability/renormalization issues before reaching any singularities, and
its current version makes no claim of singularity avoidance. Both IDG and DGR
require some early inflationary epoch to resolve the cosmological horizon
problem. Taken as a whole, IDG possesses a moderate nonlocality that does not
challenge long cherished theoretical concepts such as renormalizability. Its
elegant UV completeness arises because its form of nonlocality efficiently
suppresses high energy quantum spacetime fluctuations, and is then more
appealing to those who have conviction in renormalizability. However, so far
it does not shed much light on observational puzzles. By contrast, DGR was
explicitly designed to break the cluster decomposition principle on short
lengths so micro-causality becomes emergent rather than fundamental. [2] It
can account for some of those puzzles, yet is presently limited to being an
effective theory.

“In every journey, there are as many objectives missed as there are objectives
gained.” –Arnold Toynbee, 1960

##  11\. Acknowledgements

The author would like to thank Patrizia Vitale for bringing the work of her
research group about commutative deformations on Minkowski spacetime to his
attention. Although at first glance those integral kernel methods seemed
inappropriate to gravitational physics, in the end they indeed pointed out a
way forward.

## References

  * [1] P.G.N. de Vegvar, “Commutative deformations of general relativity: nonlocality, causality, and dark matter,” Euro. Phys. J. C 77(1) (2017) 1-26, arXiv: 1605.06011.
  * [2] Paul G.N. de Vegvar,“Who asked for causality?” in Spacetime Physics 1907-1917, Christopher Duston and Marc Holman eds., selected papers presented at the First Hermann Minkowski Meeting on the Foundations of Spacetime Physics, 15-18 May, 2017, Albena, Bulgaria, (Minkowski Institute Press, 2019).
  * [3] S. Dopplicher, K. Fredenhagen, J.E. Roberts, “The Quantum Structure of Spacetime at the Planck Scale and Quantum Fields,” Commun. Math. Phys. 172 (1995) 187; and D.V. Ahluwalia, “Quantum Measurement, Gravitation, Locality,” Phys. Lett. B 339 (1994) 301 and 303, arXiv: 9308007.
  * [4] M.R. Douglas and N.A. Nekrasov, “Noncommutative Field Theory,” Rev. Mod. Phys. 73 (2001) 977, arXiv:0106048 [hep-th].
  * [5] A.H. Chamseddine, “Deforming Einstein’s Gravity,” Phys. Lett. B 504 (2001) 33, arXiv: 009153.
  * [6] P. Aschieri, M. Dmitrijevic, F. Meyer, J. Wess,“Noncommutative Geometry and Gravity,” Class. Quant. Gravity 23 (2006) 1883, arXiv: 0510059.
  * [7] A. Schenkel, “Noncommutative Gravity and Quantum Field Theory on Noncommutative Curved Spacetimes,” Ph.D. Thesis, University of Würzburg (2012) unpublished, arXiv: 1210.1115.
  * [8] F. Lizzi and P. Vitale, “Gauge and Poincaré Invariant Regularizations and Hopf Symmetries,” Mod. Phys. Lett. A 27 (2012) 125097, arXiv:1202.1190 [hep-th].
  * [9] S. Galluccio, F. Lizzi, and P. Vitale, “Translational Invariance, Commutation Relations, and Ultraviolet/Infrared Mixing,” JHEP 0909 (2009) 054, arXiv: 0907.3540 [hep-th].
  * [10] F. Ardalan, H. Arfaei, M. Ghasemkhani, and N. Sadooghi, “Gauge Invariant Cutoff QED,” Phys. Scripta 03 (2013) 035101.
  * [11] P. Aschieri and L. Castellani, “Noncommutative Gravity Solutions,” J. Geom. Phys. 60 (2010) 375, arXiv: 0906.2774.
  * [12] M. Nakahara, Geometry, Topology and Physics, p. 186 (Institute of Physics Publishing, 2003).
  * [13] Steven Weinberg, Gravitation and Cosmology, p. 376 (John Wiley and Sons, Inc., 1972).
  * [14] M. Nakahara, op. cit., pp. 193-4.
  * [15] Lars Hörmander, The Analysis of Linear Partial Differential Operators, Vol. I, Thm. 7.1.5, p. 160 (Springer Verlag, 2003).
  * [16] Steven Weinberg, Gravitation and Cosmology, op. cit., p. 133, eqn. 6.1.5.
  * [17] M. Nakahara, Geometry, Topology and Physics, op. cit., pp 269-70.
  * [18] S. Weinberg, The Quantum Theory of Fields, Vol. II, p. 2, (Cambridge University Press, 1996).
  * [19] M. Nakahara, op. cit., Chap. 10.
  * [20] C. Rovelli, Quantum Gravity, Chap. 6 (Cambridge University Press, 2004).
  * [21] R. Gambini and J. Pullin,Loops, Knots, Gauge Theories and Quantum Gravity, (Cambridge University Press 1996).
  * [22] C. Rovelli, op. cit., chap 2.
  * [23] C. Rovelli, op. cit., pp. 145-6.
  * [24] C. Rovelli, op. cit., p. 159 Eqn. 4.68.
  * [25] A. Aropostathis, “A counter example to a nonlinear version of the Krein-Rutman theorem by R. Mahadevan,” Nonlinear Analysis, 171C (2018) 170, ArXiv: 1802.04745.
  * [26] S. Weinberg, Gravitation and Cosmology,, op. cit., p. 395.
  * [27] S. Weinberg, The Quantum Theory of Fields, Vol. I, p. 320, (Cambridge University Press, 1995).
  * [28] See ref. 1, discussion before equation 27.
  * [29] Steven Weinberg, Gravitation and Cosmology,, op. cit..
  * [30] Steven Weinberg, Gravitation and Cosmology,, op. cit., pp 62-3.
  * [31] J. Ambjørn, J. Jurkiewicz, R. Loll, Phys. Rev. Lett., 95(17) 171301 2005, arXiv: 0505113.
  * [32] L. Modesto, “Fractal Structure of Loop Quantum Gravity,” Cl. Quant. Grav., 26 242002 (2009), arXiv: 0812.2214.
  * [33] A. Addazi, P. Belli, R. Bernabei, A. Marcianò, “Testing Noncommutative Spacetimes and Violations of the Pauli Exclusion Principle with underground experiments,” Chin. Phys., C42(9) 094001 (2018), arXiv: 1712.08082.
  * [34] N.D. Birrell and P.C.W. Davies, Quantum Fields in Curved Space, p. 62 (Cambridge University Press 1982).
  * [35] S. Weinberg, The Quantum Theory of Fields, Vol. I op. cit. pp. 170-1.
  * [36] S. Weinberg, The Quantum Theory of Fields, Vol. I op. cit. pp. 173.
  * [37] Bennett et al., “Search for Lorentz and CPT Violation Effects in Muon Spin Precession,” Phys. Rev. Lett., 100(9) 091602 (2008), arXiv:0709.4670.
  * [38] CMS Collaboration,“Measurement of the mass difference between top and antitop quarks,” JHEP 06 109 (2012), arXiv: 1204.2807.
  * [39] Ch. Eisele, A. Yu Nevsky, S. Schillerv, “Laboratory Test of the Isotropy of Light Propagation to the $10^{-17}$ level,” Phys. Rev. Lett., 103(9) 090401 (2009); S. Hermann, A. Senger, K. Möhle, M. Nagel, E.V. Kovalchuk, A. Peters, “Rotating optical cavity experiment tresting Lorentz invariance at the $10^{-17}$ level,” Phys. Rev. D80 (100) 105011 (2009), arXiv: 1002.1284.
  * [40] V.W. Hughes, H.G. Robinson, V. Beltran-Lopez,“Upper Limit for the Anisotropy of Inertial Mass from Nuclear Resonance Experiments,” Phys. Rev. Lett., 4(7) 342 (1960); R.W.P Drever, “A search for anisotropy of inertial mass using a free precession technique,” Phil. Mag., 6(65) 683 (1961).
  * [41] F. Allmendinger et al., “New limit on Lorentz and CPT violating neutron spin interactions using a free precession ${}^{3}\mathrm{He}$-${}^{129}\mathrm{Xe}$ co-magnetometer,” Phys. Rev. Lett.., 112(11) 110801 (2013), arXiv: 1312.3225.
  * [42] N.V. Krasnikov, “Nonlocal gauge theories”. Theoretical and Mathematical Physics 73 (2): 1184–1190 (1987).
  * [43] E.T. Tomboulis, “Super-renormalizable gauge and gravitational theories,” arXiv: 9702146 (1997).
  * [44] T. Biswas, T. Koivisto, A. Mazumdar, “Nonlocal theories of gravity: The flat space propagator,” Phys. Rev. Lett. 108 031101 (2012), arXiv: 1302.0532.
  * [45] L. Buoninfante, G. Lambiase, Y. Miyashita, W. Takebe, and M. Yamaguchi, “Generalized ghost-free propagators in nonlocal field theories,” Phys. Rev. D 101 084019 (2020), arXiv: 2001.07830.
  * [46] J. Edholm, “UV completion of the Starobinsky model, tensor-to-scalar ratio and constraints on nonlocality,” Phys. Rev. D 95 044004 (2017), arXiv: 1611.05062.
  * [47] L. Feng, “Light bending in infinite derivative theories of gravity,” Phys. Rev. D 95 084015 (2017), arXiv: 1703.06535.
  * [48] A.O. Barvinsky, “Nonlocal action for long distance modification of gravity theory,” Phys. Lett. B572 109 (2003), arXiv: 0304229.
  * [49] M.E. Soussa and R.P. Woodward, “Nonlocal formulation of MOND,” Cl. Qntm. Grav. 20 2737 (2003), arXiv: 0302030.
  * [50] L. Modesto, “Super-renormalizable quantum gravity,” Phys. Rev. D 86 044005 (2011), arXiv: 1107.2403; L. Modesto and L. Rachwał, “Super-renormalizable and finite gravitational theories,” J. Nucl. Phys. B 10 015, arXiv: 1407.8036; L. Modesto and L. Rachwał, “Universally finite gravitational and gauge theories,” J. Nucl. Phys. B 09 006 (2015), arXiv: 1503.00261.
  * [51] S. Talaganis and A. Teimori, “Hamiltonian analysis for infinite derivative field theory and gravity,” arXiv:1701.01009.
  * [52] Y.D. Li, L. Modesto, and L. Rachwał, “ Exact solutions and spacetime singularities in nonlocal gravity,” JHEP 12 173 (2015), arXiv: 1506.08619.