# Monge Ampère gravity:
from the large deviation principle to cosmological simulations through optimal
transport

Bruno Lévy<EMAIL_ADDRESS>Inria Saclay Ile de France, Laboratoire de
Mathématiques d’Orsay, 91405, Orsay, France  Yann Brenier<EMAIL_ADDRESS>CNRS, Ecole Normale Supérieure, 45 rue d’Ulm, F-75230 PARIS cedex 05, France
Roya Mohayaee<EMAIL_ADDRESS>Sorbonne Université, CNRS, Institut
d’Astrophysique de Paris, 98bis Bld Arago, 75014 Paris, France  Rudolf Peierls
Centre for Theoretical Physics, University of Oxford, Parks Road, Oxford OX1
3PU, United Kingdom

###### Abstract

We study Monge-Ampère gravity (MAG) as an effective theory of cosmological
structure formation through optimal transport theory. MAG is based on the
Monge-Ampère equation, a nonlinear version of the Poisson equation, that
relates the Hessian determinant of the potential to the density field. We
explain how MAG emerges from a conditioned system of independent and
indistinguishable Brownian particles, through the large deviation principle,
in the continuum limit. To numerically explore this highly non-linear theory,
we develop a novel N-body simulation method based on semi-discrete optimal
transport. Our results obtained from the very first N-body simulation of
Monge-Ampère gravity with over 100 millions particles show that on large
scales, Monge-Ampère gravity is similar to the Newtonian gravity but favours
the formation of anisotropic structures such as filaments. At small scales,
MAG has a weaker clustering and is screened in high-density regions. Although
here we study the Monge-Ampère gravity as an effective rather than a
fundamental theory, our novel highly-performant optimal transport algorithm
can be used to run high-resolution simulations of a large class of modified
theories of gravity, such as Galileons, in which the equations of motion are
second-order and of Monge-Ampère type.

cosmology, optimal transport, Monge-Ampère, modified gravity, Galileons, extra
dimensions

###### pacs:

##  I introduction

The elliptic Monge-Ampère (MA) equation:

$\det(\mathrm{D}^{2}\Phi({\bf r}))=\alpha\rho({\bf r})$ (1)

is a nonlinear second-order partial differential equation, where the potential
$\Phi$ and the density $\rho$ are functions of spatial coordinates ${\bf r}$,
$\mathrm{D}^{2}$ denotes the Hessian, and $\alpha$ is a constant scalar. It
has applications in optimal transport, meteorology, optics, oceanography and
astrophysics. It exhibits some similarities with the Poisson equation,

$\Delta\Phi({\bf r})=\alpha\rho({\bf r}),$ (2)

that prescribes the sum of the eigenvalues of the Hessian
$\Delta\Phi=\mbox{trace }(D^{2}\Phi)$ whereas the Monge-Ampère equation
prescribes their product.

In one dimension, the two equations are the same. In higher dimension, the
Monge-Ampère equation can be linearized around the identity to obtain the
Poisson equation. The extra off-diagonal terms may be thought of as tidal
tensors, yielding additional anisotropic components in the associated force.
Unlike the Poisson equation, the Monge-Ampère equation is invariant under
affine symmetry which includes shear transformation. As we shall see, affine
symmetry yields anisotropic structures – e.g. ellipsoidal haloes, filaments
and sheets – more stable and abundant in Monge-Ampère gravity than in the
Newtonian gravity. Furthermore, unlike Newtonian gravity, MAG remains finite
at short distances, which imply a screening mechanism in nonlinear regimes.

The MA equation was first studied 200 years ago by Monge, and later
generalised by Ampère. In the 90’s, one of us, Brenier Brenier (1991)
established the relation with another legacy of Monge: optimal mass transport
Monge (1784), that corresponds to the gradient of the potential that solves
the Monge-Ampère equation. The existence and regularity of solutions as well
as ways of solving the Monge-Ampère equation have been active areas of
research Figalli and Glaudo (2021); Villani (2003, 2009).

In meteorology, replacing Poisson by Monge-Ampère equation was a fruitful
strategy: Eliassen’s semi-geostrophic model Eliassen (1948), provides a
successful framework to study atmospheric motion with the formation of shocks
and fronts Hoskins (1975); Cullen and Purser (1984), and works better on large
scales than the older Poisson-based quasi-geostrophic model Cullen (2006).

It has been shown by one of us, Yann Brenier Brenier (2015) that Monge-Ampère
gravity can be derived from a system of Brownian particles through the
application of large deviation principle (more on this later).

Monge-Ampère equation is also tightly related to optimal transport which has
been shown to be a powerful _inverse method_ for reconstructing the initial
condition of the Universe from the present data. The relation between the
Monge-Ampère equation, optimal transport and the least action principle
Benamou and Brenier (2000) was used to design over the years increasingly
efficient methods for early Universe reconstruction Frisch  _et al._ (2002);
Brenier  _et al._ (2003); Lévy  _et al._ (2021); von Hausegger  _et al._
(2022); Nikakhtar  _et al._ (2022, 2023), that may be thought of as variants
of Peeble’s action method Peebles (1989); Nusser and Branchini (2000).

In the present article, in contrast with the reconstruction methods mentioned
above, we study for the first time the _forward_ dynamics governed by the
Monge-Ampère equation, using high resolution numerical simulations. Here we
only consider Monge-Ampère gravity as an emergent rather than fundamental
force, at cosmological scales, where convexity holds. The validity of Monge-
Ampère gravity at smaller scales (sub Mpcs), where convexity breaks down due
to multi-streaming (Brenier  _et al._ , 2003), shall be discussed in
forthcoming works Bonnefous-Abiven  _et al._ (2024).

In the rest of this article, we first detail the definition and mathematical
properties of Monge-Ampère gravity, first based on a mathematical intuition
obtained by examining the equations (beginning of Section II) and their
relations with Optimal Transport (Section II.1). Then we re-obtain the same
equations from a physical path of reasoning, deriving the laws of motions as
the continuum limit of a microscopic model of a large number of particles
subject to the principle of large deviations, similar to a least action
principle (Section II.3). We exploit this physical reasoning to invent a new
numerical simulation algorithm (Section III.1). Finally, we run the same
N-body simulation with $512^{3}$ particles (Section IV), with Newtonian
gravity (using an adpative-mesh algorithm similar to Couchman  _et al._
(1994)) and Monge-Ampère gravity, both starting from the same $\Lambda$CDM
initial conditions, and make some qualitative observations. To our knowledge,
this is the first high-resolution simulation of highly non-linear modified
gravity, which we believe can be applied to a larger class of modified
theoriess of gravity, such as Gallileon Fairlie (2011).

##  II Monge-Ampère gravity

In this section, we show how the Monge-Ampère equation (3) relates to the
usual Poisson equation in Newtonian gravity. We then explain its relation with
Optimal Transport. We subsequently use this relation to derive Monge-Ampère
gravity from a microscopic model and to design a cosmological N-body
simulation. In Monge-Ampère gravity Yann Brenier (2004); Brenier (2011);
Ambrosio  _et al._ (2023), the gravitational force derives from a potential
$\phi$, related to the solution $\Phi$ of a Monge-Ampère equation:

$\det(\mathrm{D}^{2}\Phi({\bf r}))=\frac{\rho({\bf r})}{\bar{\rho}}\,,\\\ $
(3)

where $\bar{\rho}$ denotes the background density. Next we write

$\Phi({\bf r})=\frac{\phi({\bf r})}{4\pi{\cal G}\bar{\rho}}+\frac{|{\bf
r}|^{2}}{2}\,,$ (4)

and expand the determinant in (3) around the identity to obtain the Poisson
equation:

$\displaystyle 1+\mbox{trace }\left(\mathrm{D}^{2}\frac{\phi({\bf
r})}{4\pi{\cal G}\bar{\rho}}\right)+\mathcal{O}(\phi^{2})$ $\displaystyle=$
$\displaystyle\frac{\rho({\bf r})}{\bar{\rho}}\,,$ (5)
$\displaystyle\Delta\phi({\bf r})$ $\displaystyle\simeq$ $\displaystyle
4\pi{\cal G}(\rho({\bf r})-\bar{\rho})\,.$ (6)

This corresponds to the usual forms of the Monge-Ampère (1) and Poisson
equation (2) (using $\Phi$ as defined in eq. (3) and $\alpha=1/\bar{\rho}$).
In our Poisson equation (6), the average density $\bar{\rho}$ is subtracted
out of $\rho({\bf r})$. This is a common practice in astrophysics, since a
uniform density, $\bar{\rho}$, yields no potential, and is often neglected at
the galactic scales whereas taken into account at cosmological scales.

###  II.1 Monge-Ampère and optimal transport

Optimal Transport concerns the properties of maps $T$ that minimize the
integrated weighted quadratic transport cost,

$\inf\limits_{T}\left[\int\rho({\bf r})|{\bf r}-T({\bf r})|^{2}d{\bf
r}\right]\,,$ (7)

subject to mass conservation:

$\int_{B}\bar{\rho}d{\bf q}=\int_{T^{-1}(B)}\rho({\bf r})d{\bf r}\quad\forall
B\,,$ (8)

(for all subset $B$, there should be the same quantity of matter in $B$ and in
its pre-image through $T$). To analyze the optimization problem above, it is
easier to rewrite the constraint in terms of functions :

$\int g({\bf q})\bar{\rho}d{\bf q}=\int g(T({\bf r}))\rho({\bf r})d{\bf
r}\quad\forall g.$ (9)

Intuitively, one can see that the equation above implies (8) if one takes
$g({\bf r})$ = 1 if ${\bf r}$ is in $B$ and 0 otherwise. Enforcing a
constraint through functions $g(.)$ is a common practice in Finite Elements
Analysis (where functions $g(.)$ are referred to as _test functions_).

We shall now follow a derivation similar to the Euler-Arnold formulation of
incompressible fluids (Arnold, 1966), where pressure appears as the Lagrange
multiplier associated with incompressibility. We first replace the squared
distance in the quantity to be minimized by the dot product:

$\inf_{T}\left[\int|{\bf r}-T({\bf r})|^{2}d{\bf
r}\right]\quad\rightarrow\quad\sup_{T}\left[\int{\bf r}\cdot T({\bf r})d{\bf
r}\right]$

and we inject the constraint as a Lagrange multiplier. Then, our optimization
problem (7) subject to mass conservation (9) rewrites as a saddle-point
problem:

$\begin{array}[]{ll}\sup\limits_{T}\inf\limits_{\Psi}\scalebox{1.5}{[}&{\cal
L}(T,\Psi)=\int\rho({\bf r})T({\bf r})\cdot{\bf r}d{\bf r}\ +\\\\[5.69054pt]
&\int\bar{\rho}\Psi({\bf q})d{\bf q}\ -\int\Psi(T({\bf r}))\rho({\bf r})d{\bf
r}\ \scalebox{1.5}{]}\end{array}$ (10)

where the potential $\Psi$ is the Lagrange multiplier associated with mass
conservation (9).

From the first-order and second order optimality conditions, one deduces that:

$\displaystyle\frac{\partial{\cal L}}{\partial T}=0$ $\displaystyle\
\Rightarrow\ $ $\displaystyle{\bf r}=\nabla\Psi(T({\bf r}))$ (11)
$\displaystyle\frac{\partial^{2}{\cal L}}{\partial T^{2}}\geq 0$
$\displaystyle\ \Rightarrow\ $ $\displaystyle\Psi\mbox{ is a convex function
}$ (12)

The relation (11) above gives ${\bf r}$ as a function of $T({\bf r})$ (inverse
optimal transport map). Since $\Psi$ is convex, the (direct) optimal transport
map $T$ is given by the gradient of another convex potential $\Phi$, related
with $\Psi$ through its Legendre-Fenchel transform $\Psi^{*}$.

$\displaystyle T({\bf r})$ $\displaystyle=$ $\displaystyle\nabla\Phi({\bf
r}),\quad\mbox{where:}$ (13) $\displaystyle\Phi({\bf r})$ $\displaystyle=$
$\displaystyle\Psi^{*}({\bf r})=\inf_{{\bf q}}\left[{\bf q}\cdot{\bf
r}-\Psi({\bf q})\right].$ (14)

In thermodynamics, such convex potentials $(\Psi,\Phi)$ are said to be
_Legendre dual_. Their duality is a direct consequence of the relation between
the Lagrangian and the Hamiltonian formalisms.

Next, $T({\bf r})=\nabla\Phi$ is inserted into the mass-conservation
constraint in Eq. (9):

$\displaystyle\bar{\rho}\int g(\nabla\Phi({\bf r}))|\mathrm{D}^{2}\Phi({\bf
r})|d{\bf r}$ $\displaystyle=$ $\displaystyle\int g(\nabla\Phi({\bf
r}))\rho({\bf r})d{\bf r}$ (15) $\displaystyle\bar{\rho}\det D^{2}\Phi$
$\displaystyle=$ $\displaystyle\rho({\bf r}).$ (16)

Eq.(15) is obtained by applying the change-of-variable formula to the LHS of
Eq. (9). This yields the Monge-Ampère equation (16) by recalling that Eq. (15)
should be satisfied for all "test functions" $g(.)$.

From the above derivation, we deduce the relation between the solution $\Phi$
of the Monge-Ampère equation (3), the optimal transport map $T$ and the
gravitational potential $\phi$:

  * •

the optimal transport map is given by $T=\nabla\Phi$;

  * •

the gravitational potential $\phi$ is related to $\Phi$ by:
$\Phi=\phi/(4\pi{\cal G}\bar{\rho})+|{\bf r}|^{2}/2$, because Eq. (16) is
identical to the Monge-Ampère equation (3).

###  II.2 Symmetries of the Monge-Ampère equation

The Newtonian potential, as the solution of the Poisson equation, is rotation-
invariant, that is, if $\Phi$ is a solution of $\Delta\Phi=\rho$, then we have
$\Delta(\Phi\circ R)=\rho\circ R$ for any rotation $R$ (or rigid motion). This
explains why a spherical lump of matter (invariant by rotation) remains stable
under Newtonian gravity.

The Monge-Ampère equation shows invariance with respect to a wider class of
transforms, that preserve volumes (gobally or locally). Intuitively, the
Monge-Ampère equation can be thought of as (local) volume conservation for a
map $T$ (where the map $T$ is the gradient of a convex potential $\Phi$),
hence it is invariant with respect to a volume-conserving transform applied to
both sides of the equation. More formally, consider a solution $\Phi$ of the
Monge-Ampère equation $|D^{2}\Phi|=\rho$, and a transform ${\bf
u}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$. We shall now see under which
condition $\Phi\circ{\bf u}$ is a solution of $|D^{2}[\Phi\circ{\bf
u}]|=\rho\circ{\bf u}$. Expanding the l.h.s., one gets:

$\displaystyle|D^{2}[\Phi\circ{\bf u}]({\bf x})|$ $\displaystyle=$
$\displaystyle\left|J\left[J^{t}[\Phi\circ{\bf u}]({\bf x})\right]\right|$
(17) $\displaystyle=$ $\displaystyle\left|J[J\Phi({\bf u}({\bf x}))\times
J{\bf u}({\bf x})]^{t}\right|$ $\displaystyle=$
$\displaystyle\left|J\left[J^{t}({\bf u}({\bf x}))\times J^{t}\Phi({\bf
u}({\bf x}))\right]\right|$ $\displaystyle=$ $\displaystyle\left|\ D^{2}{\bf
u}({\bf x})\times J^{t}\Phi({\bf u}({\bf x}))+\right.$ $\displaystyle\ \
\left.J^{t}{\bf u}({\bf x})\times D^{2}\Phi({\bf u}({\bf x}))\ \ \right|,$

where $Jf$ denotes the Jacobian matrix of $f$, using the relations $\nabla
f=J^{t}f$ (gradient is a _transposed_ Jacobian), $D^{2}f=J[J^{t}f]$, the chain
rule and the product rule.

In Equation 17, $D^{2}{\bf u}$ is a $3\times 3$ tensor with coefficients that
are the second order derivatives of $u$. Hence the first term vanishes if
${\bf u}$ is an affine function, in the form ${\bf u}({\bf x})=A{\bf x}+{\bf
b}$ (or if $u$ has very small second order derivatives), then one gets:

$\displaystyle|D^{2}[\Phi\circ{\bf u}]({\bf x})|$ $\displaystyle=$
$\displaystyle\left|J^{t}{\bf u}({\bf x})\times D^{2}\Phi({\bf u}({\bf
x}))\right|$ (18) $\displaystyle=$ $\displaystyle\left|J{\bf u}({\bf
x})\right|\times\left|D^{2}\Phi({\bf u}({\bf x}))\right|.$ (19)

Hence, if ${\bf u}$ is volume-preserving, that is, with unit Jacobian, then
the term $|J{\bf u}({\bf x})|$ is equal to one, and one gets:

$|D^{2}[\Phi\circ{\bf u}]({\bf x})|=\left|D^{2}\Phi({\bf u}({\bf x}))\right|.$
(20)

To summarize, if $\Phi$ is a solution of $|D^{2}\Phi|=\rho$, consider a
function ${\bf u}$ with unit Jacobian and vanishing (or negligible) second
order derivatives, then $|D^{2}[\Phi\circ{\bf u}]|=\rho\circ{\bf u}$. In other
words, the Monge-Ampère equation is invariant w.r.t. a function that locally
behaves like a (volume-preserving) shearing. Clearly, one needs to further
analyze the implication for the Monge-Ampère gravitational potential $\phi$,
which will be done in future works. At this point, we have gained the
intuition that anisotropic structures (ellipsoids, filaments, …) are more
likely to remain stable through Monge-Ampère gravity than with Newtonian
gravity. This will be empirically illustrated in the results section.

Figure 1:  Left panel: unconditioned motion of $M$ independent Brownian
particles; Center panel: motion of independent Brownian particles conditioned
by their initial and final positions (in red and blue respectively); Right:
conditioned Brownian motion with vanishing noise, all trajectories tend to
geodesics.

###  II.3  Monge-Ampère gravity;
emergence from the large deviation principle

The derivation of Monge-Ampère gravity in the beginning of this section is
motivated by a mathematical intuition: using the expansion of a determinant
around the identity, one can construct an equation similar to Newton-Poisson
gravity with additional higher-order terms, by replacing the Laplace operator
(determinant of the trace of the Hessian) with the Monge-Ampère operator
(determinant of the Hessian), as derived in equations (3) to (6).

In this section, we detail an alternative derivation, leading to the the very
same equation, but this time motivated by a physical principle. We shall use
later this alternative physically-motivated path of reasoning to invent an
efficient numerical simulation algorithm.

Following Brenier (2015), we shall derive Monge-Ampère gravity from a
microscopic model, that shares some similarities with the relation between
optimal transport and Schrödinger bridges Léonard (2012). For the sake of
completeness, we mention that there also exists an elegant "pilot-wave"-like
derivation Benitez-Llambay (2015).

The idea resembles the least action principle: assume that a motion from a
fixed initial condition to a fixed final configuration minimizes (or
extremizes) some criterion, deduce the law of motion in-between as a
differential relation, then extrapolate it.

So we consider the time evolution of a density field $\rho({\bf r},t)$, that
accounts for a (huge) number of $M$ indistinguishable particles. The only
assumption regarding the particles is that they move randomly and do not
interact. With this only assumption, everything will spread out and converge
to a uniform distribution of particles (see the left panel of Figure 1).

Now we suppose that the density field $\rho({\bf r},t)$ was observed at time
$t=0$ and at a time $t=T$. Then we seek for the "most probable" motion from
the initial condition at time $t=0$ (red particles in Figure 1-center panel)
that accounts for the observations at time $t=T$ (blue particles in Figure
1-center panel). As can be seen, this may be thought of as a statistical
version of the least action principle, where we seek for the "most probable"
motion (instead of extremizing an action).

The initial condition at time $t=0$ is represented by a permutation
$\sigma_{0}$ of $M$ particles, that accounts for the particle’s
indistiguishability :
$\left\\{X^{0}_{1},\ldots,X^{0}_{M}\right\\}=\left\\{{\bf
q}_{\sigma_{0}(1)},\ldots,{\bf q}_{\sigma_{0}(M)}\right\\}.$ We suppose that
each particle’s trajectory is given by ${\cal X}_{i}(t)=X^{0}_{i}+\epsilon
B_{i}(t)$ where the $B_{i}(t)$’s are $M$ independent realizations of a
Brownian motion and where $\epsilon$ denotes the amount of noise.

Given a set of points $(Y_{i})_{i=1}^{M}$, the probability of the (very
unlikely, hence "large deviation") event that at time $T$, each position
$Y_{i}$ is occupied by one of the particles ${\cal X}_{i}(t)$ writes:

$\begin{array}[]{l}\mbox{Prob }\left({\cal
X}_{i}^{\epsilon}(T)\underset{\text{p}erm}{\approx}Y\right)\approx\frac{(2\pi\epsilon
T)^{-\frac{3M}{2}}}{M!}\underset{\sigma\in
S_{M}}{\sum}\exp\left[\frac{-\underset{i}{\sum}|Y_{\sigma(i)}-X_{i}^{0}|^{2}}{2\epsilon
T}\right]\par\end{array}$ (21)

where the sum is over all permutations $S_{M}$ of $[1\ldots M]$ to account for
particle’s indistinguishability.

Similarly to what happens in path integrals, among the permutations, the one
that minimizes $\underset{i}{\sum}|Y_{\sigma(i)}-X_{i}^{0}|^{2}$ has a
tremendous influence. As the noise $\epsilon$ tends to zero, the $\log$ of
this quantity (a "smoothed infimum") becomes the true infimum (Laplace lemma
or Boltzmann distribution in statistical physics):

$\begin{array}[]{l}-\underset{\epsilon\rightarrow
0}{\lim}\,\epsilon\log\mbox{Prob }\left[{\cal
X}_{i}^{\epsilon}(T)\underset{\text{p}erm}{\approx}Y\right]\approx\underset{\sigma\in
S_{N}}{\inf}\left[\frac{\underset{i}{\sum}|Y_{\sigma(i)}-X^{0}_{i}|^{2}}{2T}\right]\end{array}$
(22)

This corresponds to a discrete version of the optimal transport (7) (the
permutation $\sigma$ is the discrete version of the mass-preserving transport
$T$).

Moreover, when $\epsilon$ tends to zero, all the trajectories become geodesics
(rectilinear uniform), see Figure 1-right:

$\begin{array}[]{llcl}&{\cal
X}_{i}(t)&=&X_{i}^{0}+\frac{t}{T}(Y_{(\sigma|Y)(i)})-X_{i}^{0})\\\\[5.69054pt]
\mbox{where: }&(\sigma|Y)&=&\underset{\sigma\in S_{N}}{\mbox{Arg
inf}}\left[\frac{\sum_{i}|Y_{\sigma(i)}-X_{i}^{0}|^{2}}{2T}\right]\end{array}$
(23)

Then, one can derive a law of motion in the form of an ordinary differential
equation by using a simple property: along a uniform rectilinear trajectory
(Eq. 23), $X^{0}$ is the closest point (w.r.t the distance
$\|X-Y\|^{2}=\sum_{i}|X_{i}-Y_{i}|^{2}$) to not only
$X(T)=(Y_{(\sigma|Y)(i)})_{1}^{N})$, but also to any $X(t)$ for $t$ in
$[0,T]$. Hence, one obtains the following relation:

$t\frac{d{\cal X}_{i}(t)}{dt}={\cal X}_{i}(t)-{\bf q}_{(\sigma|{\cal
X}(t))(i)}$ (24)

Let us now see how to derive a force, that is, a second-order equation with
respect to time. First, to avoid singularity at time $t=0$, we make a change
of time variable $t=\exp(\tau)$:

$\frac{d{\cal X}_{i}(\tau)}{d\tau}={\cal X}_{i}(\tau)-{\bf q}_{(\sigma|{\cal
X}(\tau))(i)}$ (25)

Now, let us rewrite Eq. 25 with the set of all trajectories ${\cal X}=({\cal
X}_{i})_{i=1}^{N}$. The right hand side can be written as the gradient of a
function $\Psi:\mathbb{R}^{3N}\rightarrow\mathbb{R}$ of all particle
coordinates:

$\begin{array}[]{llcl}&\frac{d{\cal
X}_{i}(\tau)}{d\tau}&=&-\frac{\partial\Psi}{\partial{\bf x}_{i}}({\cal
X}_{1}(\tau),\ldots,{\cal X}_{N}(\tau))\\\\[5.69054pt] \mbox{where}&\Psi({\bf
x}_{1},..{\bf x}_{N})&=&\underset{\sigma\in
S_{N}}{\inf}\left[\frac{\sum_{i}|{\bf x}_{i}-{\bf
q}_{\sigma(i)}|^{2}}{2}\right]\end{array}$ (26)

To derive an equation that resembles Newton’s law, let us now compute the
second order derivative (acceleration). Equation 26 enjoys a very particular
property. Its right hand side remains the same when differentiating with
respect to time:

$\begin{array}[]{lcl}\frac{d^{2}{\cal
X}(\tau)}{d\tau^{2}}&=&-\mathrm{D}^{2}\left(\Psi({\cal
X}(\tau))\right)\frac{d{\cal X}(\tau)}{d\tau}\\\\[5.69054pt]
&=&\mathrm{D}^{2}\left(\Psi({\cal X}(\tau)\right)\nabla\Psi({\cal
X}(\tau))\\\\[5.69054pt] &=&\nabla\left(\frac{\|\nabla\Psi({\cal
X}(\tau))\|^{2}}{2}\right)\\\\[5.69054pt] &=&-\nabla\Psi({\cal
X}(\tau))\end{array}$ (27)

Finally, from the definition of the potential $\Psi$ (Eq. 26), it is easy to
check that the acceleration $d^{2}{\cal X}_{i}(\tau)/d\tau^{2}$ of one of the
particles is given by:

$\begin{array}[]{lcl}\frac{d^{2}{\cal
X}_{i}(\tau)}{d\tau^{2}}&=&-\frac{\partial\Psi}{\partial{\cal X}_{i}}={\cal
X}_{i}-{\bf q}_{(\sigma|{\cal X}(\tau))(i)}\\\\[5.69054pt] &=&{\cal
X}(\tau)-{\bf q}_{(\sigma|{\cal X}(\tau))(i)}\\\\[5.69054pt] &=&{\cal
X}(\tau)-\nabla\Phi({\cal X}(t))\\\\[5.69054pt] &=&-\nabla\phi({\cal
X}(\tau)).\end{array}$ (28)

In the last two lines, we consider the continuum limit, where the number $M$
of particles tends to infinity, and write the equations for the potentials
$\phi$ and $\Phi$, where $\Phi$ is the solution of the Monge-Ampère equation,
and $\phi$ is determined by $\Phi({\bf x})={\bf x}^{2}/2-\phi({\bf x})$.

To summarize the derivations above:

  1. 1.

From the principle of large deviations conditioned by both the initial
condition at $\tau=\tau_{1}$ and the observation at time $\tau=\tau_{2}$, one
deduces that between $\tau=\tau_{1}$ and $\tau=\tau_{2}$, each particle is
subject to a force;

  2. 2.

this force is equivalent to a "repulsion" by a "virtual" particle, that
corresponds to the particle’s image through optimal transport towards the
initial condition. One may be surprised that it is a "repulsive" force, and
that only a single particle is involved, but we recall that this "virtual"
particle accounts for the influence of all real particles through optimal
transport;

  3. 3.

this force derives from a potential, that is obtained as the solution of the
Monge-Ampère equation $\det(D^{2}(\Phi))=4\pi{\cal G}\rho$, similar to the
usual Poisson equation, but with additional non-linear terms (the off-diagonal
coefficients of the Hessian).

  4. 4.

the Monge-Ampère gravitational force is given by:

$\displaystyle\frac{d^{2}{\cal X}_{i}(\tau)}{d\tau^{2}}$ $\displaystyle=$
$\displaystyle{\cal X}(\tau)-{\bf q}_{(\sigma|{\cal X}(\tau))(i)}$ (29)
$\displaystyle=$ $\displaystyle{\cal X}(\tau)-\nabla\Phi({\cal X}(t))\ =\
-\nabla\phi({\cal X}(\tau))$ (30)

where the continuous version (30) uses the relation ${\bf q}=\nabla\Phi$, with
$\Phi$ the solution of the MA equation (see Equations 13 to 16).

One can gain more intuition about the repulsion from the "virtual particle"
that corresponds to the optimal-transport image. Starting from the initial
condition, and with particles that move according to $M$ independent Brownian
motions, the principle of large deviation considers the probability of a very
unlikely event: exactly one particle is observed at each position
$(Y_{i})_{i=1}^{M}$). It is known that moving along the direction of optimal
transport tends to regularize the distribution of the particles Xin  _et al._
(2016) and to maximize entropy. In our case, moving along the _opposite_
direction favors the apparition of structures, by directly leading to the
(very unlikely) event.

##  III Cosmological simulation of Monge-Ampère gravity

Figure 2: Semi-discrete simulation of Monge-Ampère gravity with the path
bundle method. To derive the force, we suppose that at a time $\tau_{0}$
matter takes the form of a uniform density field. Between times $\tau_{0}$ and
$\tau_{1}$, matter clusters into $N$ structures, symbolized as red points. In
addition, we suppose that the structures are observed at certain locations at
a later time $\tau_{2}$. From the so-conditioned large deviation principle,
one deduces the law of motion of the $N$ structures

In this section, we detail our numerical simulation algorithm for Monge-Ampère
gravity. Our goal is now to simulate a set of $N$ objects (galaxies),
represented by points, or "macro-particles" $({\bf r}_{i})$ with masses
$(m_{i})$, with initial positions ${\bf q}_{i}={\bf r}_{i}(\tau_{1})$ and
velocities ${\bf v}_{i}(\tau_{1})$.

Figure 3: Semi-discrete optimal transport as the continuum limit of discrete
optimal transport. First four columns: each macro-particle is associated with
1,9,16,25 microparticles respectively. Rightmost column: when the number of
microparticles tends to infinity, each macroparticle is associated with a
polygonal region, called a Laguerre cell, computed explicitly by our
algorithm.

Clearly one could use the discrete Monge-Ampère gravitational force (29) to
derive a N-body simulation, with $F_{i}={\bf r}_{i}-{\bf q}_{(\sigma|{\cal
X})(i)}$, where $(\sigma|{\cal X})(i)$ is the optimal permutation of the $N$
points that minimizes the (weighted) transport cost
$C(\sigma)=\sum_{i}m_{i}|{\bf r}_{i}-{\bf q}_{\sigma(i)}|^{2}$. However,
finding $(\sigma|{\cal X})(i)$ at each time-step is a very costly operation,
typically in $O(N^{2}\log(N))$ Bertsekas and Castanon (1989). This algorithmic
complexity prevents us from using discrete optimal transport in cosmological
simulations, for which we need at least $512^{3}$ particles to have sufficient
resolution to see the evolution of structures. For this reason, as explained
in the next subsection, we use instead the _continuous_ Monge-Ampère force
(30), leading to a much more efficient numerical solution mechanism that
exploits the structure of the problem (convexity and smoothness). In practice,
complexity becomes $O(N\log(N)$ Lévy  _et al._ (2021) (to be compared with
$O(N^{2}\log(N)$).

###  III.1 The path bundle method

We shall now derive the laws of motion of the $N$ macroscopic particles based
on the continuous Monge-Ampère force. To do so, we consider that each of the
$N$ macro-particles is composed of a huge number of micro-particles, with a
total number $M>>N$ of micro-particles, then we will consider the continuum
limit. For now, in the discrete setting, the mass $m_{i}$ of a macro-particle
$i$ may be thought of as the number of micro-particles in $i$.

The $M$ indistinguishable micro-particles ($M>>N$) move according to $M$
independent Brownian motions conditioned as shown in Figure 2: we suppose that
the matter that corresponds to each particle $(m_{i},{\bf r}_{i})$ at the
initial condition $\tau_{1}$ has lumped from a uniform distribution of matter,
at a "pre-initial condition" $\tau=\tau_{0}$. This uniform pre-initial
condition may be thought of as a regular grid of micro-particles.

In addition, we suppose that at time $\tau_{1}$, all the micro-particles that
clustered into a macro-particle have the same velocity. Hence, from time
$\tau_{1}$ and later, the $N$ structures remain intact and can be modeled as
$N$ macro-particles. In other words, the trajectories of the microscopic
particles that clustered into a structure form a "bundle", that merges at time
$\tau_{1}$. Since microscopic particle’s velocities match at time $\tau_{1}$,
between times $\tau_{1}$ and $\tau_{2}$, each "bundle" is tight, and appears
as a unique trajectory of a macroscopic particle.

From the large deviation principle detailed in the previous section, one
deduces that each macroscopic particle is subject to a force that results from
the combined influence of all the trajectories of the micro-particles.

Now remember the intuition that a particle is "pushed away" by its optimal-
transport image (equation 29, end of previous Section). In the present
context, at a time $\tau$ between $\tau_{1}$ and $\tau_{2}$, we consider that
a macro-particle is composed of a huge number of micro-particles. Each of
these micro-particle is "pushed away" by its optimal-transport image. The
force that the macro-particle is subject to corresponds to the sum of each of
its corresponding micro-particle’s contribution.

As shown in Figure 3, taking the continuum limit, the regular grid of micro-
particle at the pre-initial condition at $\tau=\tau_{0}$ converges to the
uniform density $\rho({\bf q},\tau_{0})=\bar{\rho}$, and the force influencing
each macro particle can be derived from the optimal transport between the
uniform density and the set of (weighted) particles $(m_{i},{\bf r}_{i})$.
Thus, the gravitational force $F_{i}$ for a given object ${\bf r}_{i}(\tau)$
writes:

$F_{i}(\tau)=-\nabla\phi(\tau)={\bf r}_{i}-\nabla\Phi({\bf
r}_{i},\tau)\quad\forall\tau\geq\tau_{1}$ (31)

where $\Phi(.,\tau)$ is the solution of the Monge-Ampère equation
$\det(D^{2}\Phi(.))=\rho(.,\tau)$.

Note that the density field $\rho(.,\tau)$ as well as the potentials $\phi$
and $\Phi$ are concentrated on the $N$ punctual masses ${\bf r}_{i}(\tau)$,
hence one needs to define the differential operators $\nabla$ and
$\det(D^{2}(.))$ in a generalized sense. In the so-called _semi-discrete_
setting, the gradient $\nabla\phi$ is replaced with the subdifferential, that
associates a region of space $V_{i}(\phi)$ to each point ${\bf r}_{i}$ Brenier
(1991); Mérigot (2011); Lévy (2015); Lévy  _et al._ (2021); von Hausegger  _et
al._ (2022); Nikakhtar  _et al._ (2022),

$\begin{array}[]{l}\nabla\phi=V_{i}(\phi)=\\\\[5.69054pt] \\{{\bf q}\ ;\ |{\bf
q}-{\bf r}_{i}|^{2}-\phi_{i}\leq|{\bf q}-{\bf r}_{j}|^{2}-\phi_{j}\quad\forall
1\leq j\leq N\\}\end{array}$ (32)

where the so-defined $V_{i}$ regions (Laguerre cells) form a partition of $V$
(Laguerre diagram), parameterized by the $N$ positions ${\bf r}_{i}(\tau)$ of
the points and by the potentials $\phi_{i}(\tau)=\phi({\bf r}_{i}(\tau))$. Put
differently, each Laguerre cell $V_{i}$ corresponds to the continuum of
optimal-transport images that "pushes" the particle $i$ (rightmost column of
Figure 3).

Integrating the influences of all the points ${\bf q}$ in $V_{i}$, one obtains
the expression of the gravitational force:

$\frac{d^{2}{\bf r}_{i}(\tau)}{d\tau^{2}}=F_{i}(\tau)={\bf r}_{i}(\tau)-{\bf
g}_{i}(\tau)$ (33)

where ${\bf g}_{i}$ denotes the barycenter of the Laguerre cell $V_{i}$. The
Laguerre cells are determined by the points $({\bf r}_{i})$ and the potentials
$(\phi_{i})$ that realize the optimal transport of the $({\bf r}_{i})$’s to
the uniform density. How to compute the potentials $(\phi_{i})$ is explained
later, in Section III.3.

As compared to the fully discrete setting considered at the beginning of this
section, each point ${\bf r}_{i}$ is repelled by the barycenter of the region
of space that corresponds to its image through optimal transport (whereas in
the discrete setting, it was repelled by a single point ${\bf
q}_{\sigma(i)}$). Not only this gives more spatial accuracy for the force, but
also it makes it possible to replace the combinatorial search for
$(\sigma|{\cal X})$ with a much more efficient smooth and convex optimization.

Figure 4: Four timesteps $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$ of a 2D
cosmological Monge-Ampère gravity simulation, using the expansion scale factor
$a$ as the time variable (3D examples shown P. 9). Each timesteps consists of
the following substeps: (1) solve the discrete Monge-Ampère equation to
compute the potentials $\phi_{i}$; (2) construct the Laguerre diagram; (3)
compute the cell barycenters; (4) compute the forces (5) update the momenta;
(6) advect the particles. The bottom row shows a zoom on the Laguerre cells of
the central zone.

###  III.2 Cosmological simulation

Based on the expression of the Monge-Ampère gravitational force above, we
shall now detail our cosmological simulation algorithm. First, we recast the
equations of motion in cosmological time and space coordinates, using the
expansion scale factor $a$ as the time variable, and co-moving coordinates
${\bf x}$ defined by ${\bf r}=a{\bf x}$:

$\frac{d{\bf x}}{da}=\frac{{\bf p}}{ma^{2}S(a)}\quad;\quad\frac{d{\bf
p}}{da}=\frac{1}{aS(a)}\left[a^{2}{\bf F_{\rm g}}+\frac{mH_{0}^{2}\Omega_{\rm
m}}{2}{\bf x}\right]$ (34)

where ${\bf F}_{g}$ denotes the gravitational force (Newtonian or Monge-
Ampère) and $S(a)=aH_{0}(\Omega_{\rm m}a^{-3}+(1-\Omega_{\rm
m}-\Omega_{\rm\Lambda})a^{-2}+\Omega_{\rm\Lambda})^{1/2},$ where ${\bf
v}=a{\bf u}$ denotes the peculiar velocities, ${\bf p}=ma{\bf v}$ momenta,
$H_{0}$ the Hubble constant, $\Omega_{\rm m}$ the matter density and
$\Omega_{\rm\Lambda}$ the cosmological constant. As we assume periodic
boundaries for our cosmological boxes, the second term in the square bracket
vanishes. For this work, our simulations assume a FLRW background.

Our simulation algorithm is very similar to a classical N-body simulation with
semi-implicit Euler time integration, the only difference is the definition of
the gravitational force ${\bf F_{\rm g}}$.

Hence, one timestep of our simulation consists of the following operations,
summarized in Figure 4:

$\begin{array}[]{ll}\ (1):&\mbox{ Solve for the potentials }(\phi_{i})\mbox{
using the points }({\bf x}_{i})\\\ \ (2):&\mbox{ Compute the Laguerre cells
}\\{V_{i}({\bf x}_{i},\phi)\\}\\\ \ (3):&\mbox{ Compute the barycenters }({\bf
g}_{i})\mbox{ of the }\\{V_{i}\\}\mbox{ cells }\\\ \ (4):&\mbox{ Compute the
forces }{\bf F}_{MA}\mbox{ using Eq. \ref{eqn:MAForce} }\\\ \ (5):&\mbox{
Update momenta: }{\bf p}\leftarrow{\bf p}+\frac{\delta a}{S(a)}a{\bf F}_{\rm
MA}\\\ \ (6):&\mbox{ Update positions: }{\bf x}\leftarrow{\bf x}+\frac{\delta
a}{ma^{2}S(a)}{\bf p}\\\ \ (7):&\mbox{ Update time: }a\leftarrow a+\delta
a\end{array}$

  * •

Step (1) solves for the potential that corresponds to the optimal transport
between the set of points $(m_{i},{\bf x}_{i})$ and the uniform density
$\bar{\rho}$, more details below;

  * •

Steps (2),(3),(4) direcly use the Laguerre diagrams computed by the KMT-Newton
algorithm. Our algorithm shares some similarities with the moving-mesh
algorithm in AREPO Springel (2010), in the sense that it also uses a dynamic
unstructured polyhedral mesh. Besides the difference that here we are
simulating a modified gravity, the main differences are that here we use the
Laguerre cells instead of Voronoi cells in AREPO and in addition our scheme is
purely Lagrangian with exact volume conservation whereas AREPO is based on an
Arbitrary Lagrangian Eulerian scheme;

  * •

Steps (5) and (6) correspond to a semi-implicit Euler scheme with the scale
factor $a$ used as a time variable, where $\delta a$ denotes the time-step;

  * •

Step (7) updates the time parameter $a$. In addition, at this step, we check
the Courant-Friedrich-Lewy (CFL) condition (Courant  _et al._ , 1928). This
condition, required to ensure the stability of the numerical simulation,
imposes an upper bound for the timestep $\delta a$. In our case, the CFL
condition imposes that during a single timestep, a particle should not travel
further away than the size of its Laguerre cell (see Gallouët and Monsaingeon
(2017) for the derivation of a CFL condition in a similar context). If the CFL
condition is not satisfied by the timestep $\delta a$, a smaller timestep
$\delta a^{\prime}$ is computed in such a way that the fastest particle meets
this condition, and several substeps $a\leftarrow a+\delta a^{\prime}$ are
computed until $a$ reaches $a+\delta a$.

Figure 5: Simulation of Monge-Ampère gravity (60 Mpc/h, $256^{3}$ particles),
and zoom on the Laguerre cells of the central region.

###  III.3 Solving the semi-discrete Monge-Ampère equation

In the N-body simulation algorithm above, step (1) requires to solve for the
gravitational potential $\phi$ in the Monge-Ampère equation (16) recalled
below:

$\det(\Phi)=\det D^{2}\left(\frac{1}{4\pi{\cal G}\bar{\rho}}\phi+\frac{{\bf
x}^{2}}{2}\right)=\frac{\rho}{\bar{\rho}},$

which is a non-linear equation. To solve it numerically, we use a Newton
method, that solves a series of linear equations that converges to the
solution.

We first give an intuition of the method using a continuous-setting derivation
(with a continuous density field $\rho$), then we shall see how to adapt it to
our semi-discrete setting (with individual point masses ${\bf
x}_{i}$,$m_{i}$).

In the continuous setting, we use the following order-1 expansion of the
Monge-Ampère operator $\det D^{2}(.)$ in function of a small update
$\delta\phi$ of the potential:

$\det D^{2}\left(\Phi+\frac{\delta\phi}{4\pi{\cal
G}\bar{\rho}}\right)\simeq\det D^{2}(\Phi)+\frac{1}{4\pi{\cal
G}\bar{\rho}}\Delta(\delta\phi).$ (35)

Plugging (35) into the Monge-Ampère equation, one gets the equations that
characterizes the update $\delta\phi$ to be applied at each Newton step:

$\det D^{2}\Phi+\frac{1}{4\pi{\cal
G}\bar{\rho}}\Delta\delta\phi\simeq\frac{\rho}{\bar{\rho}}$

or:

$\Delta\delta\phi\simeq 4\pi{\cal G}\left(\rho-\bar{\rho}\det
D^{2}\Phi\right).$ (36)

Putting everything together, one obtains the following Newton algorithm:

$\begin{array}[]{ll}\ (0):\phi\leftarrow 0\\\ \ (1):\mbox{while
}\|\rho-\bar{\rho}\det(D^{2}\Phi)\|^{*}<\epsilon\\\ \ (2):\quad\mbox{solve for
}\delta\phi\mbox{ in }\Delta(\delta\phi)=\rho-\bar{\rho}\det(D^{2}\Phi)\\\ \
(3):\quad\mbox{find descent parameter }\alpha\\\ \
(4):\quad\phi\leftarrow\phi+\alpha\delta\phi\end{array}$ (37)

We give now more details about the algorithm above:

  * •

in step (0), the potential $\phi$ is set to 0, which means $\Phi({\bf r})={\bf
r}^{2}/2$, and the transport $T=\nabla\Phi$ is initialized with the identity.
The next steps will update $\phi$, which will "morph" $T$ from the identity to
the optimal transport;

  * •

in step (1), $\|.\|^{*}$ denotes some function norm and $\epsilon$ a
convergence threshold (more on this later);

  * •

in step (2), one solves a Poisson equation to find $\delta\phi$ (one can drop
the $4\pi{\cal G}$ factor since $\delta\phi$ is used as a _direction_ of
descent). For a smooth density $\rho(.)$, in periodic space, one could use the
Fourier transform to solve the Poisson equation, as suggested in Zheligovsky
_et al._ (2010). However, as mentioned in the previous section, in our case,
the density $\rho$ is singular, and matter is concentrated on a discrete set
of points $({\bf x}_{i})_{i=1}^{N}$ with masses $(m_{i})_{i=1}^{N}$. The
potential $\phi$ becomes a vector $(\phi_{i})_{i=1}^{N}$ of values attached to
the individual points ${\bf x}_{i}$. As in the previous section, we need to
consider the differential operators acting on $\phi$ in a general sense
(subdifferentials), hence Equation (36) becomes:

$\hat{\Delta}\phi=m_{i}-\bar{\rho}|V_{i}(\phi)|$ (38)

where $|V_{i}(\phi)|$ denotes the volume of the Laguerre cell associated with
point ${\bf x}_{i}$ given the vector of potentials $(\phi_{i})$, as defined in
Equation 32. Figure 5 shows what these Laguerre cells look like in a
cosmological simulation. In Equation 38, the $N\times N$ matrix $\hat{\Delta}$
denotes the $\mathbb{P}_{1}$ Laplacian, that is, the classical Finite Element
discretization of the Laplacian projected onto a basis of piecewise linear
functions. The coefficients $c_{ij}$ of the discrete Laplacian are given by:

$\begin{array}[]{lcl}c_{ij}&=&\bar{\rho}\frac{1}{2}\frac{1}{\|{\bf x}_{i}-{\bf
x}_{j}\|}\left|V_{i}(\phi)\cap V_{j}(\phi)\right|\quad\forall i\neq
j\\\\[14.22636pt] c_{ii}&=&-\sum_{j\neq i}c_{ij},\end{array}$ (39)

where $|V_{i}(\phi)\cap V_{j}(\phi)|$ is the area of the common face shared by
the Voronoi cells $V_{i}(\phi)$ and $V_{j}(\phi)$ (or $\emptyset$ is there is
not such a face).

More formally, as detailed in Lévy (2015); Kitagawa  _et al._ (2019); Lévy and
Schwindt (2018), one can obtain the same equation by maximizing the
Kantorovich dual:

$\quad\quad K(\phi)=\bar{\rho}\sum_{i}\int_{V_{i}(\phi)}\left(|{\bf q}-{\bf
x}_{i}|^{2}-\phi_{i}\right)d{\bf q}+\sum_{i}m_{i}\phi_{i},$

then the right hand side of the discrete Poisson equation (38) is obtained
from the negated gradient of $K(.)$, and the discrete Laplacian $\hat{\Delta}$
is obtained from the Hessian matrix of $K(.)$. Interestingly, this equation
also gives a physical meaning to the stopping criterion for step (1): one can
see that the gradient of $K(.)$ – or the right hand side of the discrete
Poisson equation – correspond to the mass errors (desired mass $m_{i}$ minus
transported mass $\bar{\rho}|V_{i}(\phi)|$). We chose to stop Newton
iterations as soon as the error of the worst particle is smaller than 1% of
the particle’s mass. In other words, we use the $L_{\infty}$ norm for
$\|.\|^{*}$ in step (1).

  * •

in step(3), the descent parameter $\alpha$ is determined. We use the KMT
criterion Kitagawa  _et al._ (2019), that provably generates a sequence of
iterates that converges to the solution of the semi-discrete Monge-Ampère
equation: the descent parameter $\alpha$ is decreased until the volume of the
smallest Laguerre cell is larger than a certain threshold,
$\frac{1}{2}\left(\inf_{i}|V_{i}(0)|,\inf_{i}(m_{i})\right)$, where $V_{i}(0)$
corresponds to the Laguerre cell associated with point ${\bf x}_{i}$ with zero
potential (in other words, ${\bf x}_{i}$’s Voronoi cell).

The 3D version of the algorithm above is fully detailed in Lévy (2015), and
used in the context of early Universe reconstruction Lévy  _et al._ (2021);
von Hausegger  _et al._ (2022); Nikakhtar  _et al._ (2022, 2024). In early
Universe reconstruction, a single instance of the Monge-Ampère equation needs
to be solved, but in our case, we need to solve a Monge-Ampère equation _at
each timestep_ , and solving this highly non-linear equation is much more
costly than solving a Poisson equation involved in classical N-body
simulations: in algorithm (37), one can see that _a series of Poisson
equations_ are solved at each timestep. However, for Monge-Ampère, the CFL
condition is favorable, and authorizes much larger timesteps than with
classical N-body simulations (typically 20 times larger). Even with this
observation, running times remain much longer than a classical N-body
simulation, because each Monge-Ampère timestep takes around 40 times longer
than a classical Poisson solve or FMM time (in other words, we gain a 10x
factor and loose a 40x factor, the net factor is around 4x slower).

Figure 6: The discrete Laplacian matrix seen as a hierarchy of graphs: at
level 0, Laguerre cells are shown as white lines, and the dual graph as thick
black lines. To solve the discrete Poisson equations efficiently, the
algebraic multigrid method creates a hierarchy of coarser versions of the dual
graph (levels 1, 2 …) by "lumping" values from finer level’s nodes.

To improve computation time, we take a closer look at the discrete Poisson
equation (38) and revisit the way to solve it: in our previous works Lévy  _et
al._ (2021); von Hausegger  _et al._ (2022); Nikakhtar  _et al._ (2022), we
used a Jacobi-preconditioned conjugate gradient implemented on the GPU Buatois
_et al._ (2009). If our density fields were discretized on regular grids, it
is well known that our Poisson equations could be solved in linear time (see
Di Martino and Facciolo (2018) and references herein), using a multigrid
method. Seen in frequency space, such methods may be thought of as solving for
each frequency of the solution on a grid with a resolution well adapted to the
frequency. However, we cannot use such a multigrid method because our
discretization is supported by the cells of the Laguerre diagram, that is, an
_irregular_ mesh (white cells in Figure 6). However, there exists a class of
multigrid methods, called _algebraic multigrid methods_ (AMG), that solely
manipulate the matrix of the system, and that create a hierarchy of operators
by analyzing the sparsity pattern of the matrix, seen as a graph: the discrete
Laplacian matrix $\hat{\Delta}$ has a non-zero coefficient $c_{ij}$ each time
the Laguerre cells of ${\bf x}_{i}$ and ${\bf x}_{j}$ touch along a common
face. This defines a graph, shown in black in Figure 6. The algebraic
multigrid method creates a hierarchy of matrices, by merging some vertices of
this graph, and lumping the associated coefficients accordingly. This
hierarchy of matrices play the same role as the staggered grids in multigrid
methods, and treat each frequency with a well-adapted resolution. We used the
implementation in Demidov (2019), that uses this strategy to create a
preconditioner, and that is especially well adapted to our context, with a
large-scale Poisson system supported by an irregular grid. The overall
acceleration factor as compared to our previous implementation is between 1.5x
and 2x faster.

$\Lambda$CDM, $z=5$ Monge-Ampère, $z=5$ $\Lambda$CDM, $z=3$ Monge-Ampère,
$z=3$ $\Lambda$CDM, $z=0$ Monge-Ampère, $z=0$

Figure 7:  Comparison between simulations of $\Lambda$CDM (using an adaptive-
mesh algorithm similar to Couchman  _et al._ (1994)) and Monge-Ampère in a
cube of 300 Mpc/h, $512^{3}$ particles, z=5, 3 and 0. Projected integrated
density in a 15 Mpc/h thick slab, using a logarithmic color scale. Large-scale
similarity between the two models is striking, however MAG creates more
abundant and diffuse filaments, whereas $\Lambda$CDM creates highly-clustered
small haloes. There is weaker clustering because MAG does not diverge and is
screened at short distances. Figure 8: Power spectrum of $\Lambda$CDM and
Monge-Ampère simulations (300 Mpc/h, $512^{3}$ particles). At large scales,
$\Lambda$CDM and Monge-Ampère have the same trend (with slightly more power
for Monge-Ampère). At small scales, Monge-Ampère is screened: it has less
power than $\Lambda$CDM.

Figure 9:  "Zoom"-simulation in a cube of 60 Mpc/h, $256^{3}$ particles. Left
column: $\Lambda$CDM simulation; Right column: Monge-Ampère simulation (z=0
for both). Top row shows the projected integrated density on a 1 Mpc/h thick
slab. All images use a logarithmic color scale. In this "zoom"-simulation, one
can better observe that Monge-Ampère gravity tends to create more diffuse and
more anisotropic filaments, whereas $\Lambda$CDM creates a larger number of
small dense haloes. Moreover, voids seem emptier in Monge-Ampère gravity.

##  IV Simulation Results

We conducted numerical experiments to compare the effect of Monge-Ampère
gravity with Newtonian gravity in $\Lambda$CDM. We ran two simulations: one
with Newtonian dynamics and one with Monge-Ampère gravity, both in comoving
coordinates for an expanding Universe with periodic boundary conditions, with
$512^{3}$ particles in a box of $300$ Mpc/h with $\Lambda$CDM initial
condition and parameters from Planck collaboration Planck Collaboration and
Ade (2016). The $\Lambda$CDM simulation uses an adpative-mesh algorithm
similar to Couchman  _et al._ (1994), and the Monge-Ampère simulation uses the
algorithm described in the previous section.

Figure 7 shows a slab with a thickness of 15 Mpc/h displayed with a
logarithmic color scale, for both $\Lambda$CDM and Monge-Ampère, at redshifts
5,3,1 and 0.

The large-scale structures have globally the same shape for both simulations,
which was expected since Monge-Ampère can be considered as a small
perturbation of Poisson. At a finer scale, the Monge-Ampère simulation has
more abundant and more diffuse filaments, especially at low redshift: in
Newtonian gravity, filaments are clumpier and are transformed into multiple
spherical regions of high density. The associated power spectra are shown in
Figure 8. At small scales, Monge-Ampère has much less power than $\Lambda$CDM
(that creates dense lumps of matter). At large scales, $\Lambda$CDM and Monge-
Ampère have the same trend (with slightly more power for Monge-Ampère). This
expected since in low density regions far away from dense regions, the
nonlinearities are not important and we recover linear Poisson-like equation.
On the other hand, in high-density regions near massive bodies, the nonlinear
terms become important and, a mechanism similar to the screening mechanism in
modified theories of gravity (e.g. Vainshtein mechanism), effectively suppress
the spatial derivative of the scalar field. Hence we observe fewer dense
haloes in Monge-Ampère gravity.

| Total time (s) | timesteps | time/timestep | Newton it.
---|---|---|---|---
$\Lambda$CDM | $1.1047\times 10^{6}$ | 1797 | 10 min 14 s | N/A
MA | $1.1343\times 10^{6}$ | 95 | 40 min to 7h20 | 2 to 12
Table 1: Timings for the 300 Mpc/h simulation with $512^{3}$ particles, for
the $\Lambda$CDM and the Monge-Ampère simulations.

In terms of computation time, as summarized in Table 1, our optimized
multithreaded solver for Monge-Ampère took $1.1047\times 10^{6}$ seconds (300
hours), on a Intel Xeon Gold 5122 at 3.6 GHz, with 100 Gb of RAM. Our
$\Lambda$CDM simulation took approximately the same time, $1.1343\times
10^{6}$ seconds (monothreaded code, similar to Couchman  _et al._ (1994)).
Clearly, much faster results can be obtained for $\Lambda$CDM with recent
multithreaded codes, such as GADGET-4 Springel  _et al._ (2021) (at least 10x
faster). It is expected that our algorithm will be significantly slower than
modern $\Lambda$CDM codes, because to solve the Monge-Ampère equation, our
Newton algorithm needs to solve _multiple_ discrete Poisson equations at each
time step, between 2 to 12 in the present case (as we approach z=0, the
density fields becomes more and more irregular, and the KMT Newton algorithm
for semi-discrete optimal transport needs more and more iterations to
converge). However, the very favorable CFL condition lets us compute much
larger timesteps. This is made possible by both our unstructured
representation based on Laguerre cells, and the fact that the solutions of
Monge-Ampère equation are in general much more regular than Poisson solutions:
densities become less clustered, and velocities remain bounded. The Monge-
Ampère simulation reached $z=0$ after 94 timesteps, to be compared by the 1797
timesteps used by the $\Lambda$CDM simulation. Clearly, this very favorable
CFL coundition will not suffice to make our algorithm compete with fast
$\Lambda$CDM codes (that need solving a _single_ Poisson equation per
timestep). However, we wish to stress that computation times remain
reasonable, and our algorithm is to our knowledge the first one that can be
used in practice for high-resolution simulation of such highly non-linear
modified theories of gravity.

To see the effect of Monge-Ampère gravity at a higher resolution, we also ran
a "zoom" simulation in a cube of 60 Mpc/h, with $256^{3}$ particles (2.5x
resolution in each direction in terms of number of particles). A slab with a
thickness of 1 Mpc/h and a 3D view are both shown in Figure 9. At this scale,
one can better see the structure of the filaments, how they are fragmented by
Newtonian gravity, and how they remain continuous while more fuzzy in Monge-
Ampère gravity. In addition, one can see that the voids are emptier in Monge-
Ampère gravity (this is also visible at $z=0$ in the 300 Mpc/h simulation in
Figure 7).

Animations of both the 300 Mpc/h and 60 Mpc/h simulations are available from
http://brunolevy.github.io/videos/index.html

##  V Conclusions

In this work, we have provided the first high-resolution numerical simulation
of Monge-Ampère equation as an effective theory of gravity on large scales.
Our algorithm based on semi-discrete optimal transport theory allows a novel
high-resolution simulation of this highly nonlinear second-order differential
equation. More quantitative numerical tests of the Monge-Ampère equation shall
be presented in the forthcoming works, together with an extended algorithm for
conducting higher-resolution numerical experiments.

We shall conclude by discussing the new lines of research that follow our
results, more precisely, we ask the following question:

What is the physical significance of Monge-Ampère equation and does it have
any bearing on modified theories of gravity ? We outline three possible
scenarios in which Monge-Ampère equation can describe gravity.

In the first scenario, which we have explored in this work, MAG is not a
fundamental but an effective theory of gravity on cosmological scales. The
solar system provides one of the most stringent tests of theory of general
relativity. It is clear that MAG cannot replace Poisson equation as a
fundamental force, as the nonlinearities of the Monge-Ampère equation would
violate the solar system test Bonnefous-Abiven  _et al._ (2024). MAG yields
Poisson in the weak gravity regime but will not yield Poisson in strong
gravity regime, e.g. in the solar system. Here, we have shown that MAG can
emerge from a system of Brownian indistinguishable and independent particles
through the application of large-deviation principle. We have ran parallel
simulations with Poisson equation and Monge-Ampère equation. Comparison
between these two simulations show that MAG supresses small scale
fluctuations, and favours the formation of anisotropic structures such as
filaments. We understand this to be due to nonlinear terms that supress
clustering at small scales where MAG does not diverge. Furthermore, the Monge-
Ampère equation has a higher symmetry group, invariant through a shearing
(GL(3)), whereas the Poisson equation only has rotational invariance (SO(3)).
The affine invariance of the Monge-Ampère equation implies that the whole
range of new solutions obtained through shearing and deformation are stable.
In addition, the Laplacian, tends to homogenise the fluctuations and shift
power from large to small scales. We emphasise that, Monge-Ampère gravity is
weaker at small scales and hence yields less clustering (compare to Li  _et
al._ (2013)).

In the second scenario, MAG is a fundamental theory of gravity that describes
a scalar field which is coupled to matter with the latter satisfying the
Einstein theory of general relativity, i.e. a scalar-tensor theory. In this
scenario the Monge-Ampère equation represents the equation of motion of a
large class of Lagrangians, such as the Galileons and those related to extra-
dimensional theories, and can be considered as the fifth force (see
e.g.Fairlie (2011); Curtright  _et al._ (2012); Deffayet  _et al._ (2009);
Chow and Khoury (2009); Nicolis  _et al._ (2009); Burrage (2010); Deffayet
_et al._ (2011); Deser and Franklin (2012)). The scalar field is naturally
screened at short distances (similar to Vainshtein mechanism), i.e. at the
solar system scale, because of nonlinearies and at very large distances the
Newtonian gravity prevails. The effect of scalar-field is most prominent at
galactic scales. The physical aspects of Monge-Ampère equation as a model for
scalar fields is discussed in a companion article Bonnefous-Abiven  _et al._
(2024).

In the third scenario, MAG is a fundamental theory of dark matter whereas the
Poisson equation describes the luminous matter. Models in which the
equivalence principle between dark and luminous matter is broken, have been
studied before in different contexts (see e.g. Coc  _et al._ (2009)) and is
discussed in the companion work Bonnefous-Abiven  _et al._ (2024).

Finally, the major motivation for introducing a scalar field into theories of
gravity is to produce an alternative to dark energy and to remedy the
cosmological constant problem. Here, we have ran our simulations in a FLRW
background. The study of growth rate in MAG as an alternative fundamental
theory of gravity will be carried in the forthcoming work.

## acknowledgements

We wish to thank Mike Cullen, Gilles Esposito-Faerese, Cédric Deffayet,
Christian Léonard, André Lukas, Quentin Mérigot, Farnik Nikakhtar, Henri
Orland, Pierre Salati, Jean-Philippe Uzan and Sebastian von Hausegger for
discussions, and Denis Demidov for his help with AMGCL.

## References

  * Brenier (1991) Y. Brenier, Communications on Pure and Applied Mathematics 44, 375 (1991).
  * Monge (1784) G. Monge, Histoire de l’Académie Royale des Sciences (1781) , 666 (1784).
  * Figalli and Glaudo (2021) A. Figalli and F. Glaudo,  _An Invitation to Optimal Transport, Wasserstein Distances, and Gradient Flows_ (European Mathematical Society, 2021).
  * Villani (2003) C. Villani,  _Topics in optimal transportation_ , Graduate studies in mathematics (American Mathematical Society, Providence (R.I.), 2003).
  * Villani (2009) C. Villani,  _Optimal transport : old and new_ , Grundlehren der mathematischen Wissenschaften (Springer, Berlin, 2009).
  * Eliassen (1948) A. Eliassen, Geofys. Publ. 17 (1948).
  * Hoskins (1975) B. Hoskins, J. of the atmospheric sciences 32 (1975).
  * Cullen and Purser (1984) M. Cullen and R. Purser, J. of the Atmospheric Sciences , 1477 (1984).
  * Cullen (2006) M. Cullen,  _A Mathematical Theory of Large-Scale Atmosphere/Ocean Flow_ (Imperial College Press, 2006).
  * Brenier (2015) Y. Brenier, “A double large deviation principle for monge-ampère gravitation,” (2015), working paper or preprint.
  * Benamou and Brenier (2000) J. Benamou and Y. Brenier, Numerische Mathematik 84, 375 (2000).
  * Frisch  _et al._ (2002) U. Frisch, S. Matarrese, R. Mohayaee, and A. Sobolevskii, Nature 417 (2002).
  * Brenier  _et al._ (2003) Y. Brenier, U. Frisch, M. Hénon, G. Loeper, S. Matarrese, R. Mohayaee, and A. Sobolevskii, Mon. Not. R. Astron. Soc. 346 (2003).
  * Lévy  _et al._ (2021) B. Lévy, R. Mohayaee, and S. von Hausegger, Mon. Not. R. Astron. Soc. (2021).
  * von Hausegger  _et al._ (2022) S. von Hausegger, B. Lévy, and R. Mohayaee, Phys. Rev. Lett. 128, 201302 (2022).
  * Nikakhtar  _et al._ (2022) F. Nikakhtar, R. K. Sheth, B. Lévy, and R. Mohayaee, Phys. Rev. Lett. 129, 251101 (2022).
  * Nikakhtar  _et al._ (2023) F. Nikakhtar, N. Padmanabhan, B. Lévy, R. K. Sheth, and R. Mohayaee, Phys. Rev. D 108, 083534 (2023), arXiv:2307.03671 [astro-ph.CO] .
  * Peebles (1989) P. J. E. Peebles, Astrophys. J. Lett. 344, L53 (1989).
  * Nusser and Branchini (2000) A. Nusser and E. Branchini, Month. Not. R. Astron. Soc. 313, 587 (2000).
  * Bonnefous-Abiven  _et al._ (2024) A. Bonnefous-Abiven, Y. Brenier, and R. Mohayaee, in preparation (2024).
  * Couchman  _et al._ (1994) H. M. P. Couchman, P. A. Thomas, and F. R. Pearce, arXiv: Astrophysics (1994).
  * Fairlie (2011) D. Fairlie, Journal of Physics A Mathematical General 44, 305201 (2011), arXiv:1102.1594 [hep-th] .
  * Yann Brenier (2004) G. L. Yann Brenier, Geometric & Functional Analysis 14, 1182 (2004).
  * Brenier (2011) Y. Brenier, Confluentes Math. , 361 (2011).
  * Ambrosio  _et al._ (2023) L. Ambrosio, A. Baradat, and Y. Brenier, Analysis & PDE (2023).
  * Arnold (1966) V. Arnold, Annales de l’Institut Fourier 16, 319 (1966).
  * Léonard (2012) C. Léonard, Journal of Functional Analysis 262, 1879 (2012).
  * Benitez-Llambay (2015) A. Benitez-Llambay, “py-sphviewer: Py-sphviewer v1.0.0,” (2015).
  * Xin  _et al._ (2016) S. Xin, B. Lévy, Z. Chen, L. Chu, Y. Yu, C. Tu, and W. Wang, ACM Trans. Graph. 35, 244:1 (2016).
  * Bertsekas and Castanon (1989) D. P. Bertsekas and D. A. Castanon, Annals of Operations Research 20, 67 (1989).
  * Mérigot (2011) Q. Mérigot, Comput. Graph. Forum 30, 1583 (2011).
  * Lévy (2015) B. Lévy, ESAIM M2AN (Mathematical Modeling and Analysis) (2015).
  * Springel (2010) V. Springel, Monthly Notices of the Royal Astronomical Society 401, 791 (2010), https://academic.oup.com/mnras/article-pdf/401/2/791/3952227/mnras0401-0791.pdf .
  * Courant  _et al._ (1928) R. Courant, K. Friedrichs, and H. Lewy, Mathematische Annalen 100 (1928).
  * Gallouët and Monsaingeon (2017) T. O. Gallouët and L. Monsaingeon, SIAM Journal on Mathematical Analysis 49, 1100 (2017), https://doi.org/10.1137/16M106666X .
  * Zheligovsky  _et al._ (2010) V. Zheligovsky, O. Podvigina, and U. Frisch, Journal of Computational Physics 229, 5043 (2010).
  * Kitagawa  _et al._ (2019) J. Kitagawa, Q. Mérigot, and B. Thibert, J. Eur. Math. Soc. 21 (2019).
  * Lévy and Schwindt (2018) B. Lévy and E. L. Schwindt, Comput. Graph. 72, 135 (2018).
  * Nikakhtar  _et al._ (2024) F. Nikakhtar, R. K. Sheth, N. Padmanabhan, B. Lévy, and R. Mohayaee, arXiv e-prints (2024), arXiv:2403.11951 [astro-ph.CO] .
  * Buatois  _et al._ (2009) L. Buatois, G. Caumon, and B. Lévy, Int. J. Parallel Emergent Distributed Syst. 24, 205 (2009).
  * Di Martino and Facciolo (2018) M. Di Martino and G. Facciolo, Image Processing On Line 8, 192 (2018), https://doi.org/10.5201/ipol.2018.228.
  * Demidov (2019) D. Demidov, Lobachevskii J. of Math. 40, 535 (2019).
  * Planck Collaboration and Ade (2016) Planck Collaboration and Ade, Astron. & astrophys. 594, A13 (2016).
  * Springel  _et al._ (2021) V. Springel, R. Pakmor, O. Zier, and M. Reinecke, Monthly Notices of the Royal Astronomical Society 506, 2871 (2021), https://academic.oup.com/mnras/article-pdf/506/2/2871/39271725/stab1855.pdf .
  * Li  _et al._ (2013) B. Li, A. Barreira, C. M. Baugh, W. A. Hellwing, K. Koyama, S. Pascoli, and G.-B. Zhao, J. Cosmol. Astropart. Phys. 2013, 012 (2013), arXiv:1308.3491 [astro-ph.CO] .
  * Curtright  _et al._ (2012) T. L. Curtright, D. B. Fairlie, and H. Alshal, arXiv e-prints , arXiv:1212.6972 (2012), arXiv:1212.6972 [hep-th] .
  * Deffayet  _et al._ (2009) C. Deffayet, S. Deser, and G. Esposito-Farèse, Phys. Rev. D 80, 064015 (2009), arXiv:0906.1967 [gr-qc] .
  * Chow and Khoury (2009) N. Chow and J. Khoury, Phys. Rev. D 80, 024037 (2009), arXiv:0905.1325 [hep-th] .
  * Nicolis  _et al._ (2009) A. Nicolis, R. Rattazzi, and E. Trincherini, Phys. Rev. D 79, 064036 (2009), arXiv:0811.2197 [hep-th] .
  * Burrage (2010) C. Burrage, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 259 (2010) p. 012045.
  * Deffayet  _et al._ (2011) C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, Phys. Rev. D 84, 064039 (2011), arXiv:1103.3260 [hep-th] .
  * Deser and Franklin (2012) S. Deser and J. Franklin, Phys. Rev. D 86, 047701 (2012), arXiv:1206.3217 [gr-qc] .
  * Coc  _et al._ (2009) A. Coc, K. A. Olive, J.-P. Uzan, and E. Vangioni, Phys. Rev. D 79, 103512 (2009), arXiv:0811.1845 [astro-ph] .