# Simulation techniques for modified gravity

Claudio Llinares  Institute for Computational Cosmology,
Department of Physics, Durham University,
Durham DH1 3LE, UK
Institute of Cosmology and Gravitation,
University of Portsmouth, Dennis Sciama Building,
Portsmouth PO1 3FX, UK
Instituto de Astrofísica e Ciências do Espaço, Universidade de Lisboa,
Faculdade de Ciências, Campo Grande, PT1749-016 Lisboa, Portugal
<EMAIL_ADDRESS>
(24 April 2018; 9 October 2018; 29 October 2018)

###### Abstract

The standard paradigm of cosmology assumes General Relativity (GR) is a valid
theory for gravity at scales in which it has not been properly tested.
Developing novel tests of GR and its alternatives is crucial if we want to
give strength to the model or find departures from GR in the data. Since
alternatives to GR are usually defined through nonlinear equations, designing
new tests for these theories implies a jump in complexity and thus, a need for
refining the simulation techniques. We summarize existing techniques for
dealing with modified gravity (MG) in the context of cosmological simulations.
$N$-body codes for MG are usually based on standard gravity codes. We describe
the required extensions, classifying the models not according to their
original motivation, but by the numerical challenges that must be faced by
numericists. MG models usually give rise to elliptic equations, for which
multigrid techniques are well suited. Thus, we devote a large fraction of this
review to describing this particular technique. Contrary to other reviews on
multigrid methods, we focus on the specific techniques that are required to
solve MG equations and describe useful tricks. Finally, we describe extensions
for going beyond the static approximation and dealing with baryons.

###### keywords:

Gravitation; modified gravity; cosmology theory; dark energy; dark matter;
large-scale structure of Universe; numerical methods.

PACS Numbers: 02.60.Cb, 02.60.Lj, 02.70.-c, 04.50.Kd, 98.65.-r, 98.80.-k

##  1 Introduction

Tracking the evolution of the matter distribution in the late universe is a
multi- scale problem in both space and time. For instance, the evolution of
the largest voids in the universe have associated scales of tens of Mpc and
gigayears. On the contrary, typical scales associated with merger events in
high density regions are

kpc and megayears or event shorter. Thanks to this, studying the evolution of
the universe as a whole is a highly non-trivial task. This chapter is about
how to calculate the evolution of matter over the largest possible dynamical
range (i.e. maximizing the difference between the largest and smallest
simulated scales). While these are the most precise calculations that can be
done, they are also the most expensive and require the implementation of
complex numerical techniques.

According to the scales involved, there are four ways in which we can study
the evolution of the distribution of matter in the universe:

  * •

Background cosmology: the matter distribution is assumed to be uniform. In
this case there is no gravitational collapse and the only dynamics is
associated with the expansion rate. The resulting equations are ordinary
differential equations which can be solved analytically or with Runge–Kutta
methods.

  * •

Linear cosmology: the matter distribution is assumed to have perturbations
which are small when compared to the background density. This means that it is
possible to linearize the equations with respect to the density perturbations.
In the case of General Relativity (GR), this approximation gives rise to
ordinary differential equations for the density perturbations. By using this
approach, it is possible to obtain high accuracy predictions for the redshift
$z=0$ power spectrum for scales that correspond to frequencies $k\lesssim
0.1~{}h/$Mpc.

  * •

Nonlinear cosmology: when the matter perturbations become of order one, it is
not possible to linearize the equations with respect to the density any more.
Under these conditions, it is necessary to keep all the terms in an expansion
with respect to the density. However, as cosmological evolution is not
affected by the presence of black holes, nonlinear terms can still be
neglected in the velocities and the metric perturbations. The resulting
equations are typically partial differential equations in space and time, and
complex simulation techniques are required.

  * •

Relativistic or post-Newtonian cosmology: including relativistic species or
super-horizon scales requires the addition of higher order terms in the
Einstein’s equations. The solution of these equations can be obtained with
specialized codes that take into account post-Newtonian corrections [1, 2].

This paper will deal with the third-class of calculations in a modified
gravity (MG) context. With few exceptions, existing MG codes can simulate only
scalar tensor theories, which means that they follow only one additional
degree of freedom. This requires solving the equations of motion for this
field in addition to the standard $N$-body equations. The paper is structured
according to the class of equations that describe this field. We will describe
the methods step-by-step increasing the complexity of the equations and thus
of the required solvers. A summary of this classification, as well as the
methods that are commonly used and existing implementations are given in Table
1. An in-depth discussion of the definition of the models and, in particular,
of the recurrent concept of screening mechanism can be found in the
introductory paper of this special issue.

Dark energy models share similarities with models that are intended to explain
the dark matter component of our universe (a.k.a. MOND). For this reason,
there is also a lot to learn from the implementation of these particular
models. We summarize in Table 2 existing methods and implementations for this
family of theories.

The robustness of the methods and their implementation can be quantified by
comparing results with analytic or linear solutions. However, due to the lack
of analytic solutions in the fully nonlinear regime, different implementations
of the models are validated also by comparing codes that were written
independently by different authors. Such a comparison was presented for
instance in Refs. 2013MNRAS.436..348P and 2014A&A…562A..78L. Furthermore, a
methodical analysis of the accuracy of these codes was presented in Ref.
2015MNRAS.454.4208W. Several quantities were compared: power spectrum at
different redshifts,

Clasification of modified gravity models and their implementations. In all the
cases $L$ represents a differential operator and $f$ is a scalar function.
Type of equation Model Method References
$L(\phi,\dot{\phi},\ddot{\phi})=f(\rho,\phi)$ $L,f$ linear Dynamical DE
Modified expansion ART[6] Coupled DE Modified expansion + $G_{\mathrm{eff}}$
Refs. 2010MNRAS.403.1684B,2004PhRvD..69l3516M Non-local Modified expansion +
$G_{\mathrm{eff}}$ RAMSES [9] Vector DE Modified expansion, power spectrum and
growth history Ref. 2012MNRAS.424..699C
$L(\phi,\partial_{x}\phi,\partial^{2}_{x}\phi)=f(\rho,\phi)$ $L$ linear in
$\partial^{2}_{x}\phi$. $f(R)$ Multigrid (Jordan frame) Ref.
2008PhRvD..78l3523O ECOSMOG [12] MG-GADGET [3] Multigrid (Einstein frame) Isis
[4] Generalized chameleon Multigrid ECOSMOG [13] Symmetron Multigrid AMIGA
[14] ECOSMOG [13] Isis [4] Dilaton Multigrid ECOSMOG [13]
$L(\phi,\partial_{x}\phi,\partial^{2}_{x}\phi)=f(\rho,\phi)$ $L$ non-linear in
$\partial^{2}_{x}\phi$ DGP Multigrid with smoothing DGPM [15] Fourier with
smoothing Refs. 2009PhRvD..80f4023K, 2008PhRvD..77b4048L FFT-relaxation with
operator splitting Ref. 2009PhRvD..80j4005C Multigrid with operator splitting
ECOSMOG [19] Cubic Galileon Multigrid with operator splitting ECOSMOG [20]
Quartic Galileon Multigrid with operator splitting ECOSMOG [21]
$L(\ddot{\phi},\dot{\phi},\partial_{x}\phi,\partial_{x}^{2}\phi)=f(\rho,\dot{\phi},\phi)$
Symmetron Leap-frog Solve [22]; Isis [23] $f(R)$ Implicit ECOSMOG [24]
Disformal Leap-frog Isis [25]

velocity power spectrum, mass function, density and velocity dispersion
profile of halos and field distributions. As the base codes used for each MG
implementation were different, the comparison was made on the prediction of
the deviations from GR rather that the absolute quantities. It was found that
all the MG codes that took part in the experiment agree in these predictions
up to an accuracy of 1%. In particular, this accuracy was found for the power
spectrum of density perturbations in a range of scales that goes from below
$k=0.1$ to $k\sim 7~{}h$/Mpc, which is consistent with the scales that will be
included in the next generation of galaxy surveys.

I close this introduction by giving a warning note. Perceptive readers may
notice the absence of diagrams in the text. This is actually deliberate and
intended to force the readers to do the diagrams by themselves. Different
people will make the diagrams in different ways. This will increase the
variety of interpretations of the methods described here, which in turn, will
increase the chances of new ideas to be triggered by this modest contribution.

##  2 The $N$-Body Equations for Cold Dark Matter Simulations

A change of the gravitational theory constitutes a change of one of the
fundamental pillars on which cosmological simulations are built. Thus, before
searching for simulation methods for MG, it is important to define in a
precise way what cosmological simulations are and which equations are
involved. This is not a trivial task in the most general case. So in order to
keep the discussion focused on the changes that are required to simulate a
alternative gravitational models, and not to get lost in technicalities of
more advanced simulations, we will restrict the discussion to the simplest
type of cosmological simulations. These are simulations that deal with non-
relativistic colisionless fluid (i.e. cold dark matter). Some details of what
is necessay to simulate other kinds of fluids will be given in the last
section of this paper.

Clasification of the MOND family of modified gravity models and their
implementations.  Model Method References AQUAL Algebraic Ref.
2002MNRAS.331..909N AMIGA [27] Multigrid Ref. 1999ApJ…519..590B A…464..517T
AMIGA [30] Solve [31] Multigrid (on spherical grid) NMODY [32, 33] QUMOND
Poisson based Solve [31] Ref. 2012MNRAS.421.2598A RAyMOND [35] POR [36] Bi-
metric Multigrid Solve [31]

The simplest definition we could imagine is as follows:

######  Definition 2.1.

A cold dark matter cosmological simulation is the solution of Newtonian
equations of motion in an expanding universe for a continuous density field of
nonrelativistic collisionless matter. While this definition might be good
enough when applied to standard gravity simulations, it is important to be
more precise in the MG case. This is because there is no guarantee that the
approximations that were made to derive the Newtonian limit of GR are valid in
MG. We need to keep track of all these approximations, for which it is
necessary to derive the equations starting from basic principles. For
instance, for some specific MG models, the gravitational slip (i.e. the ratio
between time and spatial scalar perturbations of the metric) is not zero. In
these cases, new equations will have to be solved to take into account the new
degrees of freedom that may arise. The form of these new equations can be only
derived in a fully relativistic set up, (i.e. starting from the basic
principles defined by the Einstein–Hilbert action and making necessary
approximations). Thus, we re-define a cosmological simulation as follows:

######  Definition 2.2.

A cold dark matter cosmological simulation is the solution of Einstein’s
equation (or their modified versions) up to first-order in the metric
perturbations and velocities. Furthermore, the energy–momentum tensor of the
matter component is assumed to be isotropic and to represent nonrelativistic
collisionless matter. With the exception of the code presented in Ref.
2017arXiv170904703L, the integration of the equations is performed in a cubic
region of space assuming periodic boundary conditions. Given this definition,
it is now easy to derive the equations that need to be solved in the
cosmological codes. We describe the derivation for the standard gravity case
first and then describe the modifications that should be made in this
derivation when dealing with extended models.

With few exceptions, these equations are solved with $N$-body methods, which
rely on the fundamental concept of $N$-body particle. These techniques are
based on the idea of discretizing the continuous density fields and describing
them with a large set of particles. Thus, the relatively small set of
equations that define the physical model (in our case Einstein’s equations) is
substituted by a much larger set of equations that describe the trajectory of
these particles. Since there is not a one-to-one correspondence between the
$N$-body particles and the physical particles that are associated to the
continuous fields, the interactions between the $N$-body particles must be
smoothed out. See Ref. 1988csup.book…..H for a detailed discussion on the
concept of $N$-body particles and methods. During the rest of this paper we
will use the terms cosmological simulation and $N$-body simulation as
synonyms. However, the distinction between the concept of simulation and the
method that we use to run them should always be kept in mind.

###  2.1 General relativity

The aim is to solve Einstein’s equation for a cosmological nonrelativistic
collisionless fluid. In $N$-body simulations, the density of such a fluid is
described by a set of particles. So we need to re-write Einstein’s equation in
a way that its solution can be described by particles. We do this by taking
its divergence

$\nabla^{a}T_{ab}=0,$ (1)

which describes the conservation of energy. The energy–momentum tensor
associated with a set of $N$-body particles is then assumed to be the addition
of Dirac’s $\delta$ functions associated with each particle. This gives a set
of geodesics equations for each particle:

$\frac{d^{2}x_{n}^{a}}{d\tau^{2}}+\Gamma^{a}_{bc}\frac{dx_{n}^{b}}{d\tau}\frac{dx_{n}^{c}}{d\tau}=0,$
(2)

where $\Gamma^{a}_{bc}$ are the Christofell symbols, the index $n$ runs over
all the particles and $\tau$ is proper time (see for instance Exercise 1.9 in
Ref. 2004sgig.book…..C). In order to get the equations for the $N$-body codes,
we need to fix the reference frame, which can be defined by fixing the metric.
The most general perturbation of the Friedmann–Lemaître–Robertson–Walker
metric in Newtonian gauge is given by

$ds^{2}=a^{2}(t)\left[-\left(1+2\Psi\right)dt^{2}-2B_{i}dx^{i}dt+\left(1-2\Phi\right)\delta_{ij}dx^{i}dx^{j}+h_{ij}dx^{i}dx^{j}\right],$
(3)

where $t$ is conformal time and $\Phi$, $\Psi$, $B_{i}$, and $h_{ij}$ are
scalar, vector and tensor perturbations. Different codes use a different
independent variable for the integration. Common choices are supercomoving
time [40], expansion factor $a(t)$ or conformal time. The supercomoving time
has the advantage of making the hydro-dynamical equations independent of $a$,
which allows the implementation of standard algorithms defined on Minkowski
spacetime. With the exception of the works presented in Refs.
2016JCAP…07..053A and 2017JCAP…11..004A, vector and tensor perturbations are
neglected. Substituting this metric in Eq. (2) and changing the parameter of
the geodesics from proper time to conformal time gives the following equation
for the spatial coordinates of each particle:

$\ddot{x}^{i}_{n}+H(a)\dot{x}^{i}_{n}-\left(\frac{\partial\Phi}{\partial
x^{i}}\right)_{n}=0,$ (4)

where the dots represent partial derivatives with respect to conformal time
and $H(a)=\dot{a}/a$ and second-order terms in the velocities were neglected.
The treatment of the damping term in this equation (second-term in the left-
hand side) can be simplified with the following definition:

$p^{i}_{n}\equiv a\dot{x}^{i}_{n},$ (5)

which gives

$\displaystyle\dot{x}^{i}_{n}$ $\displaystyle=\frac{p^{i}_{n}}{a(t)},$ (6)
$\displaystyle\dot{p}^{i}_{n}$
$\displaystyle=a(t)\left(\frac{\partial\Phi}{\partial x^{i}}\right)_{n},$ (7)

which are the Hamilton equations that can be found in $N$-body codes. Note the
presence of expansion factors in these equations. As the universe expands,
objects become farther and farther away from each other, with a consequent
weakening of the gravitational force between them. These factors take into
account this effect. It may be surprising that the effect appears in the so-
called comoving coordinates defined by Eq. (3), in which the background
position of the objects does not change with time. Even if the objects do not
move, their interdistance (which is a covariant quantity and thus, independent
of the reference frame), still increases with time, which results in smaller
and smaller forces as time passes.111It is important not to confuse the
comoving coordinates with the comoving distance. The second is an alternative
definition of distance used by cosmologists which is defined as the distance
with respect to the distance that is calculated using the background metric.
Under this alternative definition, galaxies do not depart from each other
because of the expansion, however, the effect of the expansion is still
present in the Poisson and Hamilton equations.

In order to complete the set of equations that will be solved with the
$N$-body code, we need to find an equation for the force that appears in Eq.
(7). This will be given by the Einstein’s equations themselves222Note that I
define the Einstein’s equation as containing a cosmological constant. This may
open up the discussion about this equation representing GR or the $\Lambda$CDM
cosmological model within GR. Different authors have different opinions about
these definitions. In order to fix the terminology and to avoid getting lost
in philosophical discussions about what is GR and what is modified gravity, I
define GR as in Eq. (8). Everything else will be referred to as modified
gravity.:

$G_{ab}+g_{ab}\Lambda=T_{ab}.$ (8)

The Definition 2.2 says that we need to linearize these equations with respect
to the metric perturbations and velocities (but not with respect to the
density, which is treated nonperturbatively). 333The linearization of the
velocities mentioned here should not be confused with the linearization that
is usually made when solving the equations in the context of linear cosmology.
Once the equations are determined, not restrictions exist in the amplitude of
the velocities when solving for them with $N$-body codes. So we need to write
down the perturbed Einstein’s equations for the metric (3). In the GR case,
the $ij$ component of Einstein’s equations says that the two scalar
perturbations are equal. This is also valid in MG in the Einstein frame, but
it is not necessary valid in the Jordan frame. However, as the conformal
factor that relates both frames is very close to one, it is usually assumed
that both perturbations are equal (i.e. the present generation of simulation
codes do not include gravitational slip).

Neglecting vector and tensor perturbations gives rise to equations that depend
only on one degree of freedom. This means that all the information provided by
Einstein’s equations can be condensed into one equation. It is costumary to
obtain this information from the 00 component of Einstein’s equations. By
substituting the metric given by Eq. (3) in Eq. (8), neglecting second order
terms in $\Phi$, and substracting the zeroth-order equation (i.e. Friedmann
equation) we obtain a generalization of Poisson’s equation in an expanding
universe:

$\nabla^{2}\Phi=4\pi Ga^{2}(t)\delta\rho,$ (9)

where we have assumed the following energy–momentum tensor:

$T_{ab}=(\rho_{0}+\delta\rho)u_{a}u_{b}+(P_{0}+\delta P)(g_{ab}+u_{a}u_{b}),$
(10)

where $u^{a}$ is the 4-velocity of the fluid, $\rho_{0}$ and $P_{0}$ are
background density and pressure and $\delta\rho$ and $\delta P$ are
perturbations on these quantities (which are not assumed to be small).
Furthermore, appropriate gauge choices have to be made to remove time
derivatives of the metric perturbation. See Refs. 2011PhRvD..83l3505C,
2012PhRvD..85f3512G, 2015PhRvD..92l3517F and 2017JCAP…12..022F for a
discussion on this issue. The function $a(t)$ in Eqs. (6), (7) and (9) is
obtained as a solution of the zeroth-order Einstein’s equations (i.e. the
Friedman equations):

$\frac{\dot{a}(t)}{a(t)}=H_{0}^{2}\left[\frac{\Omega_{m}}{a(t)^{3}}+\Omega_{\Lambda}\right].$
(11)

The routines that solve Eq. (9) are some of the more complex routines in
$N$-body codes and will be the main focus of the sections that follow. Several
algorithms exist at present. For instance, particle mesh (PM) codes solve
explicitly Poisson’s equation on a grid (through multigrid or Fourier-based
methods) and estimate the force with discretization formulae for the
derivatives. The forces are then interpolated back to the particle positions
using specific smoothing kernels. The density that is required to solve the
Poisson’s equation is calculated from the particle distribution using the same
type of interpolation.444We remind the reader that in order to ensure momentum
conservation, the same kernel must be used to estimate the density on the grid
and to interpolate the force to the particles’ position. The definitions of
commonly used interpolation kernels can be found in Ref. 1988csup.book…..H.

Alternatively, in the case of direct summation codes, Poisson’s equation is
integrated once analytically and thus, the gravitational force can be obtained
by adding up individual contributions from each particle:

$(\nabla\Phi)_{n}=4\pi
G\sum_{m}\sum_{n}\frac{M_{n}M_{m}}{|x^{i}_{n}-x^{i}_{m}|^{2}}\textbf{x}_{mn},$
(12)

where the indexes $m$ and $n$ run over all the particles and $\textbf{x}_{mn}$
is a unitary vector between the particles $m$ and $n$. In practice, these very
expensive summations can be substituted by approximated summations in which
the geometry of the long range interactions is simplified. Furthermore, it is
possible to combine direct summation and PM algorithms to construct even
faster codes. Finally, there are the so-called multipole methods, which rely
on a multipolar expansion of Poisson’s equation.

Note that direct summation and multipole algorithms can be implemented only
because Poisson’s equation is linear (i.e. superposition principle holds). In
general, this property does not apply to equations that arise in MG models and
thus, direct summation, multipole or Fourier techniques are not used for
simulating MG. The exact implementation of these methods in the standard
gravity case is beyond the scope of this paper and can be found for instance
in Refs. 1998ARA&A..36..599B, 2005CSci…88.1088B, 2008SSRv..134..229D or
2011ASL…..4..204B.555Further information on how simulations work can be found
in the presentation of individual codes listed in Tables 1 and 2 and
references there in. The ultimate way of understanding how an $N$-body code
works is to write one from scratch. This is not such as big task as it seems.
State-of-the-art codes typically contain more that 100.000 lines and require
several years of development. However, minimalistic implementations for
training purposes, can consist only in 2000 to 3000 lines which can be easily
written.

To summarize, in the standard gravity case, it is possible to translate the
definition of $N$-body simulation given at the beginning of this section into
the following more practical definition.

######  Definition 2.3.

A standard gravity cold dark matter cosmological simulation is the solution of
Eqs. (6) and (7) for a large set of particles. The force that appears in these
equations can be obtained by either solving Eq. (9) on a grid (in the case of
particle-mesh codes) or using Eq. (12) (in the case of tree codes) together
with Eq. (11).

The next section reformulates this definition for the MG models that were
already implemented in codes and for which already exist nonlinear
simulations.

###  2.2 Modified gravity

Following the Definitions 2.2 and 2.3, the implementation of a MG model in an
$N$-body code has two distinctive stages. First, the $N$-body equations have
to be determined. This means, that it is necessary to find a generalization of
Eqs. (6), (7), (9) and (11). The choice of these generalized equations is not
unique. For instance, very complex theories can be mapped into simple ones.
Furthermore, in several cases it is possible to approximate complex equations
with simple ones for which simpler methods can be applied without loss of
accuracy. The decisions taken at this stage can make the difference between
being able to simulate a model with minimal effort and not simulate it at all.
Once the equations have been fixed, it is possible to implement them in
$N$-body codes. In many cases the resulting equations largely differ from the
type of equations that can be found in GR. However, these equations may exist
in different fields of physics. Multidisciplinary collaborations and
communication between different communities are crucial when it comes to the
implementation stage.

####  2.2.1 Scalar field theories

Almost all the theories that have been simulated can be defined by adding a
scalar field to the Einstein–Hilbert action. There is almost infinite freedom
in the way this can be done. However, all this diversity can be condensed into
four different properties of the action. So the most general action that has
been simulated is

$\displaystyle S=$
$\displaystyle\int\sqrt{-g}Rd^{4}x+\int\sqrt{-g}\left[f(\partial\phi,\partial^{2}\phi)+V(\phi)\right]d^{4}x+$
$\displaystyle\int\sqrt{-\tilde{g}}\mathcal{L}_{M}(m(\phi),\tilde{g}_{ab}(g_{ab},\phi))d^{4}x=S_{G}+S_{KG}+S_{M},$
(13)

where $S_{G}$ is the action that describes the geometry, $S_{KG}$ the one that
gives the Klein–Gordon equation for an uncoupled field and $S_{M}$ is the
matter action. The four characteristics that define different classes of
theories are as follows:

  * •

$f(\partial\phi,\partial^{2}\phi)$: a free function of first and second
covariant derivatives of the field $\phi$ (i.e. the kinetic term), which
includes contractions between derivatives.

  * •

$V(\phi)$: a potential in which scalar field oscillates or roles.

  * •

$m(\phi)$: a coupling between the scalar field and the matter fields.

  * •

$\tilde{g}_{ab}(g_{ab},\phi)$: a relation between the Jordan and Einstein
frame metrics.

Different choices of these four ingredients give different families of
theories. The classification that is relevant for numericists is the one that
highlights the numerical complexity of the resulting equations. In this sense,
the exact form of the free functions does not matter, but only its functional
dependence. The complexity varies from very simple models in which the scalar
field is uncoupled and only affects the background evolution of the Universe
up to theories that introduce nonlinear couplings with the matter fields which
can lead to fully nonlinear partial differential equations in space and time.
Table 2.2.1 summarizes this classification with increasing complexity of the
equations and methods. Note that this classification is inspired in both the
original definition of the models as well as the numerical techniques that

Classification of scalar field models that has been simulated.  Class Kinetic
term Potential Coupling Jordan metric Model 1 $(\partial\phi)^{2}$ $V(\phi)$ 1
$g_{ab}$ Ratta&Peebels SUGRA 2 $(\partial\phi)^{2}$ $V(\phi)$ $m(\phi)$
(nonuniversal) $g_{ab}$ Coupled DE 3 $(\partial\phi)^{2}$ $V(\phi)$ 1
$A(\phi)g_{ab}$ Chameleon Symmetron Dilaton 4 $(\partial\phi)^{2}$ $V(\phi)$ 1
$A(\phi)g_{ab}+B(\phi)\phi_{,a}\phi_{,b}$ Disformal 5
$(\partial\phi)^{2}+\mathrm{higher~{}order}$ $V(\phi)$ 1 $g_{ab}$ Galileon

have been applied. However, the classification has some flexibility: theories
from a given class may be mapped into theories that belong to different
classes (see for instance Ref. 2013PhRvD..87h3010Z for a detailed discussion
on the relation between classes 4 and 5 in Table 2.2.1). Implementation
techniques that are applied to each class follow.

Class 1. Unperturbed fields: The simplest imaginable scalar fields are the
ones that have a standard kinetic term (i.e. $f=(\partial\phi)^{2}$) and are
uncoupled. As the field is not affected by metric nor matter perturbations, it
is expected that it will not have spatial perturbations and thus, it will only
evolve in the background. The only effects that this will produce in the
matter distribution are through the energy of the field, which will affect the
expansion of the Universe. Two modifications have to be made to the codes to
include this class of models: background expansion tables (which are needed
because of the dependence of Eqs. (6), (7) and (9) on $a$) and the power
spectrum that is required for generating the initial conditions (i.e. the
growth factor). Implementation of these models has no overhead with respect to
GR simulations and thus it is possible to simulate these models with the same
resolution that can be reached in state-of-the-art standard gravity
simulations.

Class 2. Fields with explicit coupling to dark matter: The simplest way to add
a coupling to matter is to couple only the dark matter component and leave the
baryons uncoupled. This will allow the model to pass solar system tests and at
the same time allow the use of very simple equations which do not include
screening mechanisms. The consequence that such coupling has in the $N$-body
equations is the addition of an effective gravitational constant which is a
function of the field plus an additional friction term.

Class 3. Fields with conformal coupling: Giving a field dependence to the
Jordan frame metric provides an extra term in the conservation of energy which
acts as a fifth force in the geodesics equation (Eqs. (6) and (7)). In
addition to this, the equation of motion for $\phi$ becomes nonlinear.
Different choices of $V(\phi)$ and $A(\phi)$ give rise to different nonlinear
terms which in turn give rise to different screening mechanisms. Thanks to
these screening mechanisms, these fields can be coupled to the baryons and at
the same time pass solar system tests. The standard kinetic term gives rise to
a hyperbolic Klein–Gordon equation which can be solved with the methods
described in Sec. 4. However, it is customary to assume the quasi-static
approximation (i.e. to neglect the time derivatives of the field). This will
give rise to elliptic equations and enable us to implement the multigrid
methods described in Sec. 3.

Class 4. Fields with disformal coupling: The addition of derivatives of the
field in the definition of the Jordan frame metric introduces a coupling with
the time derivatives of the field in both the Klein–Gordon equation as well as
in the geodesics. This means that the equations contain terms such as
$\dot{\phi}^{2}\rho$ and $\dot{\phi}^{2}\nabla\phi$. These terms can only be
treated with the nonstatic solvers described in Sec. 4.

Class 5. Fields with non-standard kinetic terms: The addition of higher-order
terms in the kinetic part of the scalar field action can give rise to fully
nonlinear equations of motion for the field (i.e. the equations are nonlinear
in the second-derivatives of the field). These equations were included in the
codes in the quasi-static limit using the multigrid techniques presented in
Sec. 3. The full nonlinearity makes the additional techniques presented in
Sec. 3.5.3 mandatory. Note that these models do not include coupling to matter
(both explicitly or through the definition of the Jordan metric). This means
that the geodesics equations are unchanged and the effects of the scalar field
on the matter distribution are through the energy of the scalar field which
will appear as an additional term in the right-hand side of the Poisson’s
equation (Eq. (9)).

Theories such as $f(R)$ are not included in Table 2.2.1 because, even if it is
possible to map them into scalar-tensor theories, their definition have
different motivation (i.e. they are not defined in the Einstein frame by
adding a scalar field, but by adding a function of the Ricci scalar to the
Einstein–Hilbert action in the Jordan frame). Details on these other theories,
are given below in independent sections.

####  2.2.2 Nonlocal gravity

A particular representation of the model can be defined by adding a term
$m^{2}R\Box^{-2}R$ to the Einstein–Hilbert action[50, 51, 52] and was already
implemented in the Ramses code[9]. It can be shown that this nontrivial
nonlocal action, can be recast into a local action which includes two scalar
fields. The resulting equations of motion for these fields are linear and
thus, it is possible to use the simple methods that are applicable to class 2
described in Sec. 2.2.1.

####  2.2.3 The $f(R)$ family of theories

The action that defines this family of theories is usually given in the Jordan
frame:

$S_{f(R)}=\int\frac{\sqrt{-\tilde{g}}}{16\pi
G}\left[R+f(R)\right]d^{4}x+S_{m}(\tilde{g}_{ab},\psi),$ (14)

where $f$ is a free function and $\psi$ are matter fields. The only model that
has been simulated to date was proposed in Ref. 2007PhRvD..76f4004H and it is
defined with the following function

$f(R)=-m^{2}\frac{c_{1}\left(\frac{R}{m^{2}}\right)^{n}}{c_{2}\left(\frac{R}{m^{2}}\right)^{n}+1}.$
(15)

The linearized Einstein’s equation includes two scalar degrees of freedom,
which means that in addition to a modified Poisson’s equation, the model
requires an extra equation to track this extra field. These two degrees of
freedom can be parametrized

as the usual gravitational potential $\Phi$ (i.e. the scalar perturbation of
the metric defined by Eq. (3)) and the derivative of the free function itself
$f_{R}\equiv\frac{df}{dR}$. After minimizing the action (14) and assuming the
quasi-static limit, we obtain the following linearized equations for these two
degrees of freedom:

$\displaystyle\nabla^{2}\Phi$ $\displaystyle=\frac{16\pi
G}{3}\delta\rho-\frac{1}{6}\delta R(f_{R})$ (16)
$\displaystyle\nabla^{2}f_{R}$ $\displaystyle=-\frac{8\pi
G}{3c^{2}}\delta\rho+\frac{1}{3c^{2}}\delta R(f_{R}),$ (17)

where $\delta R(f_{R})$ is the perturbation of the Ricci scalar with respect
to the background cosmological value. Its functional dependence with the field
$f_{R}$ is nonlinear and given by

$\delta
R(f_{R})=\bar{R}(a)\left(\sqrt{\frac{\bar{f}_{R}(a)}{f_{R}}}-1\right),$ (18)

where $\bar{R}(a)$ and $\bar{f}_{R}(a)$ are background quantities, the free
parameters of the model were fixed such that the expansion is consistent with
observations (i.e. such that it gives an effective cosmological constant) and
$n=1$. Under these conditions, the model has only one free parameter $f_{R0}$,
which sets the normalization of the background field $\bar{f}_{R}(a)$. The set
of equations that are necessary to describe nonlinear cosmology is completed
with the geodesics equations. As these theories do not include additional
fields besides the metric, the geodesics are unchanged and thus, the
trajectory of the particles is determined by $\nabla\Phi$. By combining Eqs.
(16) and (17) it is possible to see that the MG contribution to the
acceleration is given by $c^{2}\nabla f_{R}/2$.

Equations (16) and (17) substitute Eq. (9) in the Definition 2.3. These two
equations must be solved together, however, the second equation is independent
of the metric perturbation $\Phi$ and thus, can be solved independently to
obtain $f_{R}$ in every time step of the simulation. After this is done, the
field $f_{R}$ can be introduced in Eq. (16) to obtain the metric perturbation
and thus, the force that is required to move the particles.

A different implementation can be constructed by taking into account that the
theory can be mapped into a theory that includes a chameleon field[54]:

$S_{f(R),E}=\int\frac{\sqrt{-g}}{16\pi
G}\left[R+(\delta\phi)^{2}+V(\phi)\right]d^{4}x+S_{m}(\tilde{g}_{ab},\psi),$
(19)

where $g_{ab}$ and $\tilde{g}_{ab}$ are Einstein and Jordan frame metrics
respectively. This mapping can be done by defining a scalar field $\phi$ such
that

$\exp\left(-\frac{2\beta\phi}{M_{\mathrm{P}}}\right)=f^{\prime}(R)\equiv
f_{R},$ (20)

where $M_{\mathrm{P}}$ is the Planck mass and the coupling constant is
$\beta=1/\sqrt{6}$. The relation between the Einstein and Jordan frame metrics
is

$\tilde{g}_{ab}=\exp\left(-\frac{2\beta\phi}{M_{PL}}\right)g_{ab}.$ (21)

In the new frame (i.e. the Einstein frame where $f(R)=0$), the equations for
the two degrees of freedom (Eqs. (16) and (17)) are substituted with the
following:

$\displaystyle\nabla^{2}\Phi$ $\displaystyle=4\pi G\delta\rho-\rho_{f_{R}}$
(22) $\displaystyle\nabla^{2}f_{R}$
$\displaystyle=-\frac{1}{a}\Omega_{m}H_{0}^{2}\left(\eta-1\right)+a^{2}\Omega_{m}H_{0}^{2}\left[\left(1+4\frac{\Omega_{\Lambda}}{\Omega_{m}}\right)\left(\frac{f_{R0}}{f_{R}}\right)^{\frac{1}{n+1}}-\left(a^{-3}+4\frac{\Omega_{\Lambda}}{\Omega_{m}}\right)\right],$
(23)

where $\rho_{f_{R}}=\dot{f}_{R}^{2}+|\nabla
f_{R}|^{2}+V_{\mathrm{eff}}(f_{R})$ is the energy of the scalar field (which
is assumed to be small) and $\eta$ is the matter density in terms of the
background density.

As the part of the action (19) that defines the geometry is the same as in GR,
the MG effects on the mater distribution appear as a fifth force in the
geodesics equations:

$\frac{d^{2}{\bf x}}{d\tau^{2}}+\nabla\tilde{\Phi}+\frac{1}{2}\nabla f_{R}=0,$
(24)

where $\tau$ and $\tilde{\Phi}$ are super-comoving time and potential, which
are the variables that are used in the only code that includes this
representation of the model [4].

Once a frame is chosen, we need to find a way of solving the field equations.
In both cases, the field $f_{R}$ satisfies the equation of a chameleon field.
This is a quasi-linear equation, which can be solved with the nonlinear
multigrid methods described in Sec. 3. In order to ensure that the field is
positive and to increase the stability of the solvers, it is customary to
solve the equations for a transformed field $u$, which is defined by

$f_{R}=-\frac{A}{a(t)^{b}}e^{u},$ (25)

where the normalization $A$ and the dependence $b$ with the expansion factor
was chosen differently by different authors. While this change of variables
provides robustness to the solvers, it forces them to evaluate an exponential
function in every point of the grid and every iteration step. This is a very
expensive function to calculate and thus the implementation of this change
reduces the performance of the solvers. A workaround to this problem, which is
based on a different change of variables that enables the implementation of an
explicit solver, is described in Ref. 2017JCAP…02..050B.

The two formalisms described above are mathematically equivalent. However, the
numerical errors associated with the solution of these two set of equations
may be different and thus, there is no guarantee that the end result of
simulations ran in the Jordan or Einstein frame will be the same. This has
been tested in Ref. 2015MNRAS.454.4208W, were three different implementations
of the model were compared (two of them written in the Jordan frame and one in
the Einstein frame). The differences between the predictions provided by these
codes are below 1%, which is enough for the accuracy required by present and
near future surveys.

####  2.2.4 Models with higher number of dimensions (DGP)

Models defined in higher dimensional space-times, such as the DGP model, [56]
can be mapped into the scalar tensor theories described in Sec. 2.2.1, which
include the Vainshtein screening mechanism (class 5 in Table 2.2.1) [15, 57,
19]. Once this is done, it is possible to apply the same techniques applied to
fully nonlinear equations (see Sec. 3 and in particular Sec. 3.5.3).

####  2.2.5 Vector fields

The only vector-tensor theory that has been simulated is defined by the
following action [10]:

$S=\int\sqrt{-g}\left[-\frac{R}{16\pi
G}-\frac{1}{4}F_{ab}F^{ab}-\frac{1}{2}\left(\nabla_{a}A^{a}\right)^{2}+R_{ab}A^{a}A^{b}\right]d^{4}x,$
(26)

where666Note that a different sign convention was used by these authors.

$F_{ab}=\partial_{a}A_{b}-\partial_{b}A_{a}.$ (27)

The equations of motion for the four extra degrees of freedom are

$\square A_{a}+R_{ab}A^{b}=0.$ (28)

As the part of the action that defines the geometry is the same as in GR,
Einstein’s equations are unchanged, with the exception of the addition of the
energy–momentum tensor of the vector field in the source. This extra term does
not add numerical complexity. However, the equation of motion for the vector
(Eq. (28)) may be extremely difficult to solve in the most general case.

In the only existing simulations[10], the perturbations of the vector field
were not tracked and only the effects on the background expansion and the
linear cosmology were taken into account. This means that the equations Eqs.
(6), (7) and (9) were kept as described in Definition 2.3. Only the solution
of Eq. (11) was changed while running the simulation. In addition to this, the
initial conditions were modified to take into account the presence of the
vector field at high redshift. These modifications apply to both the initial
power spectrum as well as the growth rate. In other words, these authors made
appropriate approximations which map the theory into the class 1 defined in
Table 2.2.1 of Sec. 2.2.1.

####  2.2.6 The MOND case

From the point of view of late time cosmology, it is possible to classify MG
theories into two large families which are motivated as alternative solutions
to the two main open problems in cosmology: dark matter and dark energy. While
this special issue focuses in the second class of theories, the equations
associated with both families of theories share similarities. Thus, there is a
lot to learn from numerical implementations of solutions of the MOND
equations, which are representative of alternatives to dark matter.

While there are several relativistic extensions, the original MOND idea is
defined in a non-relativistic set up [58]. The theory substitutes the usual
Newtonian potential with a modified potential, which is solution of

$\nabla\cdot\left[\mu\left(\frac{|\nabla\Phi_{M}|}{a_{0}}\right)\nabla\Phi_{M}\right]=4\pi
G\delta\rho,$ (29)

where $\mu$ is a free function and $a_{0}$ a free parameter which can be
fixed, for instance, by requiring the theory to explain rotation curves of
galaxies in the absence of dark matter. A classical action for this equation
can be found in Ref. 1984ApJ…286….7B. The equation is quasi-linear and thus
can be solved with the same multigrid methods that will be described in Sec.
3.

An alternative method commonly used to solve the MOND equation consists in
equating the left-hand side of Eq. (29) with the left-hand side of the
standard Poisson’s equation. Integrating once gives

$\mu\left(\frac{|\nabla\Phi_{M}|}{a_{0}}\right)\nabla\Phi_{M}+\nabla\times\textbf{k}=\nabla\Phi,$
(30)

where the curl term k is an integration constant (in the sense that its
divergence is equal to zero). It can be shown that this term is exactly zero
for specific symmetries and that behaves at least as $r^{-3}$ for nonsymmetric
configurations[59]. By assuming that the term is second-order, it is possible
to obtain an algebraic equation for the spatial gradient of the modified
potential, which is straightforward to solve once the Newtonian potential is
calculated with the usual methods that are implemented in GR codes. The method
is exact for 1D or 2D codes respecting certain symmetries.

Note that while this method is approximate, it conserves the two main features
of the MOND theories: forces decaying as $1/r$ at large distances and a
screening mechanism based on the gradient of the field.777In fact, the MOND
community does not use the expression “screening mechanism”. Instead, they
refer to the “Newtonian limit” of the theory. However, these two terms make
reference to the same thing: the fact that these theories include mechanisms
that can hide the MG effects when needed. In the MOND case, the Newtonian
limit (where the MG effects disappear) is obtained when the gradient of the
MOND potential is large. The validity of the approximation for the highly
nonsymmetric configurations associated with the cosmic web was tested in Ref.
2008MNRAS.391.1778L (unpublished results associated with this paper can be
found also in

Ref. 2011PhDT…….253L). Among other things, it has been found that the curl
term has the effect of boosting structure formation at intermediate scales.
Since this result was obtained in cosmological models that do not include dark
matter (the dark matter effect is given by MOND), extrapolation of these
results to the dark energy case must be made with care. See Sec. 3.3 for a
description of how this approximate solution can be used as an initial guess
for iterative solvers of Eq. (29).

An alternative method is based on a different realization of the MOND theory
(QU-MOND)[60]. First simulations with this model were presented in Ref.
2011PhDT…….253L. See also other implementations presented in Refs.
2012MNRAS.421.2598A, 2015MNRAS.446.1060C and 2015CaJPh..93..232L. The
gravitational potential of this theory is described by the following set of
equations:

$\displaystyle\nabla^{2}\Phi_{N}$ $\displaystyle=4\pi G\delta\rho$ (31)
$\displaystyle\nabla^{2}\Phi_{M}$
$\displaystyle=\nabla\cdot\left[\nu\left(\frac{|\nabla\Phi_{N}|}{a_{0}}\right)\nabla\Phi_{N}\right],$
(32)

where $\Phi_{N}$ is the Newtonian potential which is solution of the standard
Poisson’s equation, $\Phi_{M}$ is the MOND potential which will affect the
trajectory of the particles and the $\nu$ free function is not to be confused
with the $\mu$ interpolating function required in the AQUAL model. Note that
this representation of the MOND idea is straightforward to implement: once the
Newtonian potential was obtained by solving Eq. (31) with usual methods, it
can be discretized on a grid and used as source of the Poisson’s equation
associated to the MOND potential (32), which is also linear and can be solved
with the same methods used for solving Eq. (31).

A different representation of the MOND idea can be derived from a bi-metric
Lagrangian and was presented in Ref. 2009PhRvD..80l3536M. This particular
theory can be simulated by combining solutions of standard Poisson’s solvers
and nonlinear solvers of Eq. (29). A particular implementation can be found in
Ref. 2011PhDT…….253L.

####  2.2.7 Code units

The reference frame adopted in simulations is usually defined by the metric
given by Eq. (3) (i.e. comoving spatial coordinates with some particular
choice of time). On top of this, it is still necessary to choose units. It is
customary to describe the space and time coordinates in terms of the size of
the box and the age of the universe ($1/H_{0}$). The unit of mass is usually
such that the mass of the particles is order one or that the mean density of
the universe is exactly one.

In the MG case, it is also necessary to chose units for the extra degrees of
freedom. These are usually defined in terms of the background value of the
field, its vacuum value or the Planck mass. Rescaling the fields with a given
power of the expansion factor may also help in removing some terms from the
Klein–Gordon or geodesics equations.

Note that these units are not natural units and thus, the speed of light and
Planck constant are not equal to one in these simulations. They must be
included explicitly with the correct units or in the normalization of the
fields. This happens for instance when taking into account time derivatives in
the equations for the extra degree of freedom.

####  2.2.8 A word on performance

Modified gravity simulations usually need to track additional degrees of
freedom, so more calculations are likely to be required. This will result in
an overhead with respect to GR cosmological simulations, which in some cases
may make MG simulations infeasible. This overhead will naturally depend on the
methods that are implemented, which in turn depends on the class of models we
are dealing with. The simplest models that can be simulated (i.e. those that
only modify the background expansion or are condensed in a rescaling of the
gravitational constant) have a straightforward implementation, which consists
in substituting the tables that have the information about the background
cosmology or the units associated to the force. This kind of simulations has
zero overhead and thus, can be run with the same resolution used for standard
gravity simulations. The same applies to nonlocal gravity (Sec. 2.2.2) when
appropriate approximations are made.

The situation is different for models which give rise to nonlinear elliptic
equations (second- and third-classes of models presented in Table 1). This
kind of equations are typically solved with iterative methods. The convergence
rate strongly depends on the specifics of both the models and the solvers. See
for instance Ref. 2017JCAP…02..050B for an example of how different
integration variables can have a dramatic change in the convergence rate of
the solvers. The overhead associated to these kind of simulations typically
lie between a factor of 2 to 10 or 20.

Finally, the most complex models we can deal with give rise to hyperbolic
equations (fourth-class of models in Table 1). The only exact solvers that
exist at present are conditionally stable, which means that they require a
Courant–Friedrich–Levy condition on the time step. The number of time steps
associated to the MG solver increases by typically two to three orders of
magnitude with respect to the standard gravity part of the code.
Implementation of implicit solvers for this family of models will largely
increase the competitiveness of these simulations in terms of speed.

##  3 Multigrid Methods for Nonlinear Partial Differential Equations

Most of the models summarized in the previous section are described by
nonlinear hyperbolic equations (i.e. nonlinear wave equations). This is
different to what happens in the GR case, where an appropriate choice of the
gauge can hide the time derivatives of the metric perturbations. In order to
keep the spirit of standard GR simulations (in which the metric perturbations
are not evolved in time starting from some initial conditions, but calculated
at every time step) it is usually assumed that the terms that include time
derivatives in the MG equations are sub-dominant or important only in the
background. This is the so-called quasi-static limit of the models, in which
the additional degrees of freedom are described by elliptic equations. The
nonlinearity of these equations prevents us from applying techniques that rely
on the superposition principle (such as Fourier-based solvers or direct
summation of forces). Instead, codes take advantage of multigrid techniques,
which are the subject of this section.

To fix the notation, let us assume that we are interested in solving an
equation that has the following form:

$\displaystyle L(\partial^{2}\phi,\partial\phi,\phi)$
$\displaystyle=S(\rho,\phi)$ $\displaystyle\mathrm{in}~{}\Omega$ (33)
$\displaystyle\phi(x)$ $\displaystyle=f(x)$
$\displaystyle\mathrm{in}~{}\partial\Omega$ (34)

where $L$ is a differential operator (not necessarily linear), $S$ is a
nonlinear scalar function of the density and the scalar field itself, $f$ is a
function that defines the boundary conditions, $\Omega$ is the domain where
the equation is defined and $\partial\Omega$ the boundary of $\Omega$. We are
interested in finding a numerical solution of this equation on a grid that
covers the simulation domain. The method consists in discretizing the equation
on the grid and using an iterative method to obtain improved solutions
starting from an initial guess. In this section we will summarize different
aspects of the iteration method as well as multi-level techniques for
accelerating the convergence rate.

The multigrid techniques are described in great detail in Refs. Brandt77,
Wesseling92, Trottenberg and brandt2011multigrid. Furthermore, a lot of
information and tricks needed to solve MG equations can be found in the papers
associated with specific implementations of the models summarized in Tables 1
and 2. Finally, as these methods are widely used in several other contexts,
such as fluid dynamics, geophysics, medicine, etc., information may be found
also in nonastrophysical journals.

###  3.1 Gauss–Seidel iterations

Before setting up an iteration scheme, it is necessary to discretize Eq. (33)
on the grid. Such discretization takes the following form:

$L^{l}[\phi^{l}]=S^{l}[\rho^{l},\phi^{l}].$ (35)

The index $l$ represents the grid in which the problem is described. The
fields $\rho^{l}$ and $\phi^{l}$ are the density and the required solution on
the grid. Finally, $L^{l}$ is a discretization of the differential operator
$L$. Equation (35) represents a very large set of nonlinear algebraic
equations, whose number is given by the number of nodes in the grid. Our aim
is to approximate the solution of the differential equation (33) by the
solution of the algebraic equation (35).

The discretization of the equation can be obtained by discretizing each
derivative of the differential operator, writing them as linear combinations
of values of the field at the neighbouring nodes. For instance, second-order
formulas for the second derivatives can be obtained using a 19-points stencil:

$\displaystyle\partial_{x}^{2}\varphi_{i,j,k}$ $\displaystyle=$
$\displaystyle\frac{1}{h^{2}}\left(\varphi_{i+1,j,k}+\varphi_{i-1,j,k}-2\varphi_{i,j,k}\right)$
(36) $\displaystyle\partial_{x}\partial_{y}\varphi_{i,j,k}$ $\displaystyle=$
$\displaystyle\frac{1}{4h^{2}}\left(\varphi_{i+1,j+1,k}-\varphi_{i+1,j-1,k}-\varphi_{i-1,j+1,k}+\varphi_{i-1,j-1,k}\right),$
(37)

where the indexes ${i,j,k}$ corresponds to the indexes of each grid in each
Cartesian direction.

Equation (35) is solved iteratively with a Gauss–Seidel scheme. In every
iteration step of this scheme, the solution is updated in every element of the
grid by finding the root of the following function

$T^{l}[\phi^{l}]\equiv L^{l}[\varphi^{l}]-S^{l}.$ (38)

In the linear case, it is possible to find an analytic expression for the
updated value $\bar{\phi}^{l}$, which is a linear combination of the values of
the field in the neighbouring nodes. This gives rise to the so-called explicit
solvers. For instance, the 1D standard Poisson’s equation, can be discretized
as follows:

$\frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{h^{2}}=\rho_{i}.$ (39)

In an explicit solver the solution in a given iteration step can be updated
for instance in the following manner:

$\tilde{\phi}_{i}=\frac{1}{2}\left[\phi_{i-1}+\phi_{i+1}-h^{2}\rho_{i}\right].$
(40)

Exact details will depend on the type of discretization and exact algorithm
used, but in general, each iteration step consists in updating the solution in
every grid point $i$ with a linear combination of the values that surround $i$
and that correspond to the previous step.

In the nonlinear case, it is customary to use implicit solvers, which consist
in approximating the root of $T^{l}$ by applying one step of a Newton–Raphson
algorithm:

$\bar{\phi}^{l}=\phi^{l}-\frac{T^{l}}{\frac{\partial
T^{l}}{\partial\phi^{l}}}.$ (41)

Reference 2017JCAP…02..050B has shown that implicit solvers are not the only
route to deal with nonlinear problems. It is in fact possible to obtain
explicit algorithms by making appropriate changes of variables.

In order to fully specify a Gauss–Seidel scheme, we need to fix the order in
which the iterations are made in each sweep through the grid. For instance, in
the lexicographic ordering, the update of the field in the cell $(i,j,k)$ is
followed by a calculation in the cell $(i+1,j,k)$ and so on. A better
convergence rate and paralelization properties can be obtained by altering
this scheme. In the most popular scheme, each complete iteration step consists
in two sweeps in which the update of a given cell $(i,j,k)$ is followed by an
update on the cell $(i+2,j,k)$. The two sweeps start in different cells:
$(0,0,0)$ and $(1,0,0)$ and thus, the scheme resembles a two colour chess
board. Similar schemes exists which consist in four or even eight colours.

###  3.2 Multigrid acceleration

The quality of a given solution $\phi^{l}$ is usually determined by studying
the amplitude of the residuals, which are defined as

$\epsilon^{l}=L^{l}[\phi^{l}]-S^{l}[\rho^{l},\phi^{l}].$ (42)

The smaller the residual, the better the solution will be. An associated
quantity, which can be used to gauge the quality of the given solver, is the
so-called convergence rate, which can be defined as the ratio of residuals at
two consecutive iteration steps. It can be shown through analysis of the
convergence rate of individual Fourier modes [62, 63, 64, 65] that the
Gauss–Seidel algorithm is very efficient in smoothing out the residuals, but
has a very slow convergence rate when the solution includes modes whose
wavelength is much larger than the spatial resolution of the grid. This is due
to the local nature of the method: in each iteration step, information travels
from each cell to its closest neighbours only and thus a very large number of
iterations is required to transfer information across the domain.

It is possible to accelerate the convergence rate for these long modes by
combining solutions obtained on coarser grids and coarser grids.888These
different grid levels should not be confused with the different grid levels
that exist in codes that include adaptive mesh refinements (AMR). In that
case, different grids are used to increase the resolution locally, while in
the multigrid case described here, resolution is actually downgraded below the
target resolution where we want to obtain the solution. The typical way of
arranging these grids is to define a set of grids whose resolutions are one
half, one quarter and so on of the target resolution. The coarsest grid of
this set usually contains two nodes per dimension, but this is not a strict
rule. The coarsest grid can even be changed while the iterations progress. In
order to prove convergence properties of the method, it is enough to work with
two grids. The algorithm is the following: a number of iterations is performed
on the finest grid $l$ (i.e. the grid in which we want to obtain the
solution). When the convergence rate slows down (or a fixed number of
iterations has been completed), iterations are made to approximate the error
of this intermediate solution. As most of the error will come from the large
modes that did not converge yet, this is done on the next coarse grid $l-1$.
So the new equation for which iterations are made on the coarse grid is

$L^{l-1}[\delta\varphi^{l-1}]=R(\epsilon^{l}),$ (43)

where $R$ is a restriction operator which transfers fields from the fine to
the coarse grid. Once a given number of iterations is made on the coarse grid,
the additional

information that was collected related to the large modes is added to the
solution on the fine grid:

$\bar{\phi}^{l}=\phi^{l}+P(\delta\phi^{l-1}),$ (44)

The transfer between the coarse and fine grids is made with a prolongation
operator $P$. See for instance Refs. Brandt77,Wesseling92,Trottenberg or
brandt2011multigrid for a detailed description of convergence properties of
this method.

Application of the method to real problems requires repeating this process on
more than one grid level. By doing this, it is possible to reconstruct the
solution on the finest grid by combining solutions of coarser and coarser
grids which treat each set of Fourier modes separately. The decision about
changing the grid level can be made by taking into account the behaviour of
the residuals and convergence rate or by using fixed schemes. Popular schemes
are the V cycle (in which the iterations are made starting from the finest
grid, all the way down to the coarsest grid and then back to the finest grid)
or the W cycle in which jumps between coarse levels are more often that
between fine levels. Well behaved solvers typically reduce the residual by one
order of magnitude after each V cycle.

In the nonlinear case, it is not possible to add up solutions obtained in
different grids and thus, the iterations in the coarse levels are not made to
approximate the error $\delta\phi^{l-1}$, but the solution $\phi^{l-1}$
itself. So Eq. (43) is substituted by

$L^{l-1}[\varphi^{l-1}]=-R(\epsilon^{l}(\varphi^{l},S^{l}))+\epsilon^{l-1}(R(\varphi^{l}),R(S^{l}))$
(45)

The correction then is made in the following way:

$\bar{\varphi}^{l}=\varphi^{l}+P(\varphi^{l-1}-R(\varphi^{l})).$ (46)

This extension to the method is known with the name of full approximation
storage (FAS).

###  3.3 Initial guess for the iterative solvers

Making appropriate choices for the initial guess with which the iterations
start can be crucial for reaching convergence in a reasonable number of
iterations. In general, this choice depends on the particulars of the equation
to be solved. The decision is usually made by trial and error. The initial
guess for obtaining solutions on the domain grid (i.e. the grid that covers
the entire simulation box) can be obtained in the following ways:

  * •

Background value of the field.

  * •

Minimum of the effective potential.

  * •

Extrapolation in time from previous time step.

  * •

Solution of auxiliary PDEs. This method has been applied to the MOND equation
(29) in Ref. 2011PhDT…….253L (see Appendix A in this reference). An auxiliary
PDE can be constructed by taking the divergence of the field $\nabla\Phi_{M}$,
which can be obtained after neglecting the curl term $\nabla\times\textbf{k}$
in Eq. (30). This is a linear Poisson’s equation and thus, can be solved with
Fourier-based algorithms which are both efficient and easy to implement. When
using this solution as the initial guess, the nonlinear multigrid solver that
is required to solve Eq. (29) does not need to make iterations to obtain the
main features of the solution, but only needs to recover the effects
associated with the curl field, which are small.

The initial guess for obtaining solutions in the refined levels is usually
determined by interpolating values from the closest coarse grid. This
naturally applies only to AMR codes.

###  3.4 Convergence criteria

An important decision that needs to be made when developing a multigrid solver
is related to the criterion that is used to stop the iterations. Criteria that
are too permissive may give inaccurate solutions, while criteria that are too
strict may force the solver to make many more iterations than required to
understand the underlying physics and may become extremely inefficient
(especially in the context of $N$-body simulations, where the equations must
be solved several hundreds or even thousands of times in a single simulation).
When designing a convergence criterion, there are three different solutions to
take into consideration:

  * •

$\phi^{l}_{n}$ = approximate solution of the discretized equation (35) on the
grid $l$ which was obtained after $n$ iterations.

  * •

$\phi^{l}$ = exact solution of the discretized equation (35) on the grid $l$.

  * •

$\phi$ = exact solution of the PDE (33).

The following are important facts associated with these solutions:

  * •

The multigrid method converges after infinite number of iterations to the
correct solution of the discretized equation:

$\phi^{l}_{\infty}=\phi^{l}.$ (47)

  * •

The fact that the convergence rate $\eta$ in the grid $l$ is close to one
(which means that the solution does not evolve when making iterations) does
not imply that the iterations have converged (it only implies that the
iterations for the modes associated with the grid $l$ have converged):

$\left(\eta\equiv\frac{\phi^{l}_{n}}{\phi^{l}_{n+1}}\sim
1\right)\nRightarrow\left(\frac{\phi^{l}_{n}}{\phi^{l}_{\infty}}\sim
1\right).$ (48)

  * •

The solution of the discretized equation is not the solution of the
differential equation evaluated at the nodes of the grid:

$\phi^{l}_{\infty}\neq\phi.$ (49)

This is because the discretization of the equation introduces an error, which
depends on the order of the discretization formulas used and on the resolution
reached by the grid. A commonly used estimation of the truncation error
$\tau^{l}$ is999The truncation error can be estimated by comparing the result
of applying the discretized differential operator in two consecutive grids $l$
and $l-1$ to the analytic solution of the equation:
$\tau^{l}=L^{l-1}\left[\phi\right]-L^{l}\left[\phi\right]$. In practice, the
analytic solution is substituted by the only solution that we have, which is
the fine grid solution $\phi^{l}$. In order to be able to apply the coarse
grid differential operator $L^{l-1}$ to the fine solution $\phi^{l}$, we need
to restrict it by applying the restriction operator $R$. This gives rise to
the first term of the right-hand side of Eq. (50). Similarly, once we applied
the differential operator in both grids, we need to restrict the fine one to
be able to make the comparison (which gives rise to the second-term of the
right-hand side of this equation).

$\tau^{l}=L^{l-1}\left[R(\phi^{l})\right]-R(L^{l}\left[\phi^{l}\right]).$ (50)

  * •

The residual is not a dimensionless quantity, which means that it is possible
to make it arbitrarily small by changing the units:

$[\epsilon^{l}]\neq 1$ (51)

  * •

The residual is not necessarily uniform on the grid and may indeed have large
variations across the domain:

$\epsilon^{l}(\phi^{l}_{n})=\epsilon^{l}(\phi^{l}_{n},x)$ (52)

This is because the convergence rate is not uniform. More complex features in
the solution in specific regions will require more iterations.

Convergence criteria are usually based on amplitude of the residual, so the
first thing we need to do is to fix is a norm for the residual. Common choices
are norm 2 or $\infty$:

$\displaystyle||.||^{l}_{2}$
$\displaystyle=\Sigma_{i,j,k}\left(\epsilon_{i,j,k}^{l}\right)^{2}$ (53)
$\displaystyle||.||^{l}_{\infty}$ $\displaystyle=\max|\epsilon_{i,j,k}^{l}|,$
(54)

where the indices $(i,j,k)$ run over all the grid points in the 3D box.

A few commonly adopted criteria follow:

  * •

Comparison with initial guess:

$\frac{||\epsilon_{n}||}{||\epsilon_{0}||}<a.$ (55)

This criterion has the problem that if the initial solution is already good,
then the solver will make more iterations than necessary. It is possible to
also have the opposite case, in which the initial solution is too bad. In this
case, reducing the residual by a given fraction of the initial residual, may
not be good enough. The technique is implemented in several solvers with
standard or MG. However it must be used with care.

  * •

Comparison with a fixed constant:

$||\epsilon_{n}||<a.$ (56)

This criterion has the problem that $a$ must be chosen according to the units
used in the code.

  * •

Comparison with truncation error:

$\frac{||\epsilon_{n}||}{||\tau_{n}||}<a.$ (57)

In this case, the iterations will stop once the error associated with the
discretization of the equations dominates.

In all these criteria, $a$ is a small free parameter.

###  3.5 Extended multigrid techniques

Details of the solver strongly depend on the properties of the equation that
we are dealing with. Convergence properties depend on both the differential
operator and the source of the equation. In some cases it may be necessary to
apply additional tricks to reach convergence in a reasonable number of
iterations. We describe a few examples.

####  3.5.1 Selective iterations

In situations in which the complexity of the solution strongly depends on the
position in the domain (such as for instance systems in which the density is
almost uniform in a large part of the domain and has strong variations in a
small region) the solvers do not usually converge uniformly in the whole
domain. In some regions the solution will quickly approach the desired
solution, while in others, it may require a much larger number of iterations
to reduce the local residual. In this cases, it is possible to implement
selective iterations, which means that some local criterion is established to
stop the iterations. Thus, the iterations will be made only in the region in
which the solution has not converged and thus the overall computational effort
will be reduced.

####  3.5.2 Continuation method

The continuation method can be applied to nonlinear equations for which a
particular parameter defines a linear limit [62]. In these cases, it may be
convenient to solve the linear equation in the first place and obtain
intermediate solutions changing the values of the parameter such that the
solution evolves from the linear case towards the fully nonlinear situation in
a continuous way.

####  3.5.3 Dealing with fully nonlinear equations: Operator splitting

Models that include Vainstein screening mechanism such as DGP or Galileon give
rise to fully nonlinear partial differential equations for the additional
degrees of freedom. The equations are still second-order in space and time,
but there is a nonlinear dependence with the second spatial derivatives, which
represent a major challenge for multigrid solvers. Naive multigrid
implementations such as those described in Secs. 3.1 and 3.2 have poor
convergence properties in high density regions, which make simulations
impractical or simply impossible to run. We describe here two different
strategies that are commonly used to deal with this problem. A comparison
between different Vainstein solvers can be found in Ref. 2015MNRAS.454.4208W.

The first attempts to solve these equations implemented workarounds to this
problem which consisted in using a spherically symmetric approximation (which
gives equations that can be solved with a Fourier-based algorithm)[57, 17] or
by smoothing out the density distribution before sending it to the multigrid
solver[15, 57, 16]. The second approach can improve the convergence rate and
make the fully nonlinear simulations feasible. However, it comes with a price
in accuracy.

An alternative method, the operator splitting technique, was proposed in Ref.
2009PhRvD..80j4005C and implemented afterwards in the code ECOSMOG[19, 21].
Let us fix notation by assuming that we are interested in finding solutions of
the following equation:

$\nabla^{2}\phi-A\left[\left(\nabla^{2}\phi\right)^{2}-(\partial_{i}\partial_{j}\phi)(\partial^{i}\partial^{j}\phi)\right]=B\delta\rho,$
(58)

which has the same form of the DGP equation of motion and where $A$ and $B$
are model dependent constants and $\partial_{i}$ are partial derivatives with
respect to spatial coordinates. In order to define the Gauss–Seidel
iterations, we will need to solve a discretized version of the equation for
the field $\phi_{l,m,n}$ in each node $(l,m,n)$. If we can ensure that the
discretization of the term
$(\partial_{i}\partial_{j}\phi)(\partial^{i}\partial^{j}\phi)$ in the node
$(l,m,n)$ is independent of $\phi_{l,m,n}$, then we can see the equation as an
equation for $(\nabla^{2}\phi)_{l,m,n}$ and solve it analytically for
$\phi_{l,m,n}$ once we obtained $(\nabla^{2}\phi)_{l,m,n}$. This can be done
by decomposing the term $\partial_{i}\partial_{j}\phi$ into trace and
traceless parts:

$\partial_{i}\partial_{j}\phi=T_{ij}+S_{ij},$ (59)

where

$\displaystyle T_{ij}$ $\displaystyle=\frac{\nabla^{2}\phi}{3}I$ (60)
$\displaystyle S_{ij}$
$\displaystyle=\partial_{i}\partial_{j}\phi-\frac{\nabla^{2}\phi}{3}I,$ (61)

where $I$ is the identity operator. Substitution of these definitions in
$(\partial_{i}\partial_{j}\phi)(\partial^{i}\partial^{j}\phi)$ gives

$(\partial_{i}\partial_{j}\phi)(\partial^{i}\partial^{j}\phi)=\left(\frac{\nabla^{2}\phi}{3}\right)^{2}+\left[2S_{ij}T^{ij}+S_{ij}S^{ij}\right],$
(62)

which in turn gives the following alternative version of the equation of
motion (58):

$\nabla^{2}\phi-A\left\\{\frac{2}{3}\left(\nabla^{2}\phi\right)^{2}-\left[2S_{ij}T^{ij}+S_{ij}S^{ij}\right]\right\\}=B\delta\rho.$
(63)

The term inside the square brackets in this expression fulfils the property
that we are interested in: the usual discretization that can be obtained in
the cell $(l,m,n)$ by applying Eq. (37) is independent of the value of the
solution $\phi_{l,m,n}$ in that cell.

It has been found that this way of arranging the terms largely improves the
convergence rate of multigrid solvers. However, it does not solve all the
issues with this model since the solution of Eq. (63) for $\nabla^{2}\phi$ may
give rise to square roots of negative numbers in under-dense regions, where
density perturbations are negative. This does not occur in the DGP case, but
in the quartic Galileon[21, 66] model. See Ref. 2015PhRvD..92f4005W for a
detailed discussion on this issue.

###  3.6 Summary

This section summarizes the decisions that need to be taken when developing
multigrid codes. This list of decisions should not be considered as a simple
summary, but as a checklist of properties of the solver that can be revised
when things do not go as expected. So if the solver that the reader is
developing has poor convergence rate, it may be worth checking the following
properties of his or her implementation:

  * •

Parametrization of the theory: Can we find an alternative representation of
the theory which simplifies the equations?

  * •

Integration variables (independent variables as well as fields).

  * •

Stencil and discretization formula for the PDE inside the domain.

  * •

Stencil and discretization formula for the PDE at the boundaries.

  * •

Initial guess for the iterations.

  * •

Type of iterations (explicit, implicit, over-relaxation).

  * •

Sweeping strategy (2, 4, 8 colors?).

  * •

Number of multigrid levels.

  * •

Number of iterations per multigrid level (which may be different from level to
level).

  * •

Multigrid scheme ($V$, $W$, adaptive, etc.).

  * •

Restriction and prolongation operators.

  * •

FAS versus no FAS scheme.

  * •

Convergence criterion.

  * •

Extended techniques: selective iterations, continuation methods, operator
splitting, etc.

Taking these decisions constitutes a form of art, whose outcome depends on the
particular details of the equations to be solved. What works for one equation
may not work for others.

Once a newly implemented solver can pass several tests and is accepted as
being accurate enough, it may be a good idea to optimize it. A few ways in
which this can be done are discussed in Refs. 2015JCAP…12..059B and
2017JCAP…02..050B. The first of these optimization methods consists in
reducing the spatial resolution when solving the Klein–Gordon equation for the
additional degree of freedom. The impact of resolution in these solutions was
also discussed in Ref. 2014PhRvD..89h4023L.

In cases in which the code needs to solve both, the Poisson’s equation for the
metric perturbation and the Klein–Gordon equation for the scalar field, it is
possible to increase the speed of the simulations by solving both equations in
parallel. As the GR and MG solvers work with very similar routines, it is
possible to share the loops that have to be made through the grid to calculate
different quantities. This method was implemented in the code Isis [4] for the
symmetron model and it was found to give an improvement in speed of about 20%.

##  4 Nonstatic Solvers for Nonlinear Partial Differential Equations

The dynamics of MG degrees of freedom is usually described by nonlinear
hyperbolic equations (please, see Refs. weinberger1995first or
agarwal_equations, your favourite book on partial differential equations or
the notes of your undergrad courses if in need of refreshing classification of
equations). A common practice is to assume the so-called quasi-static
approximation, which consists of assuming that the time derivatives of these
fields are small or important only in the background. This approximation
usually gives rise to nonlinear elliptic equations, which can be solved with
the methods described in Sec. 3. This section is about going beyond this
approximation.

###  4.1 The quasi-static approximation

Terms that include time derivatives of the fields are not exclusive to the
field equations associated with MG. For instance, in the GR case, it is
necessary to deal with time derivatives when deriving the usual Poisson’s
equation from basic principles. The effects associated with these terms were
analysed in Refs. 2011PhRvD..83l3505C,2012PhRvD..85f3512G,2015PhRvD..92l3517F
and 2017JCAP…12..022F.

The validity of the quasi-static approximation in a MG context was discussed
for first time in Ref. 2008PhRvD..78l3523O. In that early work, the time
derivatives were not take into account when running the simulations, but
estimated afterwards by post-processing the simulation data (i.e. by applying
appropriate discretization formulas which combine static solutions that were
obtained at different time steps). The results show that the terms that
include time derivatives in the $f(R)$ equations are sub-dominant when
compared with the spatial derivatives. A similar comparison was made by myself
using symmetron simulations that were presented in Ref. 2013PhRvL.110p1101L. I
found that, while it is true that globally the time derivatives are sub-
dominant, they may be larger than the spatial derivatives when restricting the
analysis to under-dense regions (i.e. voids), where the Laplacian of the field
is very small.

These results can only be confirmed by running simulations beyond the static
limit and comparing the data with the output of static codes, which do not
rely of the quasi-static approximation. This has been done in the linear
regime for the $f(R)$ model[71] and in the nonlinear regime for two different
models: symmetron and $f(R)$.[22, 23, 24, 5]. The approximation was found to
be good enough for predicting the matter and velocity power spectrum. For
instance, in the case of the symmetron model, a difference of 0.2% was found
between static and non-static power spectra in a range of scales from $k$=0.05
to $k$=10 $h$/Mpc[23]. Differences in the mass function were studied in Ref.
2015MNRAS.454.4208W and also found to be well below 1% for masses between
$3\times 10^{12}$ and $10^{14}$. In the case of $f(R)$ model, differences
between static and nonstatic results are within the same ranges[24]. However,
it is worth keeping in mind that these nonlinear codes contain
simplifications. Furthermore, the impact of these additional terms on
different observables is still unknown.

It is worth noticing that there are models, such as disformal gravity, in
which the time derivatives of the field appear not only in the kinetic part of
the equations of motion, but also in the coupling to matter.[25]. This
provides a very rich phenomenology which does not exist in models that are
conformally coupled. Neglecting time derivatives when running simulations with
this kind of model may prevent us from predicting the unique signature of new
physics that we are looking for. So nonstatic simulations are crucial for
testing these models.

A problem that arises when running nonstatic simulations is related to the
determination of initial conditions for the additional degrees of freedom.
Models that include screening mechanisms are expected to give rise to small
perturbations at the start of the simulations. However, the exact shape and
amplitude of the spectrum of these perturbations may still be important. This
issue was discussed in Ref. 2016A&A…585A..37H, where it was shown that the end
result of the simulations is independent to the amplitude that is chosen for
the initial perturbations of the field. In this particular model (disformally
coupled symmetron), this is not only a consequence of the screening, but also
of the fact that there is a phase transition at a redshift of about one. The
scalar field receives a large kick when this happens, which forces it to loose
memory of any small perturbations that may have existed at early times. Note
that this result may not be valid for models that do not include screening of
the field at high redshift. In this case, initial conditions will have to be
generated more carefully, for instance, as a realization of a power spectrum
generated with linear theory.

###  4.2 Solving hyperbolic equations in $N$-body simulations

The first solvers for hyperbolic equations in the context of MG $N$-body
simulations were applied to the symmetron model and are based on the leap-frog
algorithm[22, 23]. The method consists in converting the second order
Klein–Gordon equation into a first order system of equations with the
following change:

$q=a^{b}\dot{\phi},$ (64)

where $b$ is chosen such that commonly existing terms cancel. The resulting
equation has the following form:

$\displaystyle\dot{\phi}$ $\displaystyle=\frac{q}{a^{b}}$ (65)
$\displaystyle\dot{q}$ $\displaystyle=S(\rho,\phi),$ (66)

where $S$ is a model dependent function. Note that these equations are
schematic and that details may depend on the model. These equations are
similar to Hamilton equations for a set of particles with the difference that
instead of evolving the position of the particles, they evolve the amplitude
of the field at every point of space. So by discretizing the field on a grid,
it is possible to apply the same algorithm that is commonly used to integrate
Hamilton’s equations in $N$-body codes. The method was applied to uncoupled
scalar fields within GR for instance in Refs. 1989PhRvD..40.1002W and
1989ApJ…347..590P. Implementation to coupled scalar fields can be found in
Refs. 2013PhRvL.110p1101L, 2014PhRvD..89h4023L and 2016A&A…585A..37H. As in
the model studied in these papers the scalar field undergoes very fast
oscillations, the simulation has to deal with two different time scales (one
associated with normal gravity and a much shorter one associated with the
scalar field). This problem has been solved by having two different time steps
in the simulation: a main time step as in a usual GR code and a much smaller
one for modified gravity. As the time step that is employed to move the
particles is the one associated with GR, these simulations still have the
limitation that the particles can not see the oscillations of the scalar field
and so the impact of the scalar waves in the matter distribution is still an
open question. The same method has been implemented in 1D in spherical
coordinates[67, 74] or 2D cartesian grid[75].

The leap-frog algorithm has the disadvantage that it is conditionally stable,
which means that it requires a Courant–Friedrich–Levy condition on the time
step: For the method to work, the time step has to be such that no information
is propagated in more than one spatial cell element within the step[76, 77].
While this condition does not affect the quality of the results, it forces the
simulation to work with extremely small time steps, which makes them very
expensive. A work around to this problem consists of introducing numerical
viscosity which suppress small scale oscillations of the solution and can be
done by means of implicit methods.

A different implementation of nonstatic solvers was presented in Ref.
2015JCAP…02..034B and applied to the $f(R)$ model. In this case, the time
derivatives of the field were estimated with a down-wind discretization:

$\displaystyle\left(\frac{\partial\phi}{\partial t}\right)_{n}$
$\displaystyle=\frac{\phi_{n}-\phi_{n-1}}{t_{n}-t_{n-1}},$ (67)
$\displaystyle\left(\frac{\partial^{2}\phi}{\partial t^{2}}\right)_{n}$
$\displaystyle=\frac{({d\phi}/{dt})_{n}-(d\phi/dt)_{n-1}}{t_{n}-t_{n-1}},$
(68)

which involves information from a given time step $n$ and the previous one
$n-1$. Once this is done, it is possible to translate the terms that depend on
these derivatives to the right-hand side of the Klein–Gordon equation:

$\nabla^{2}\phi_{n+1}=f(\phi_{n+1},\rho_{n+1})+g(\dot{\phi}_{n},\ddot{\phi}_{n}),$
(69)

where $f$ and $g$ are model dependent functions. This way of discretizing the
Klein–Gordon equation gives rise to an elliptic equation for the scalar field
at the time step $n+1$ and, which can be solved with static solvers described
in Sec. 3. Detailed discussion of up-wind and down-wind discretizations can be
found for instance in Refs. TORO:60926 and Press:2007:NRE:1403886.

Additional information about integration of nonstatic equations in gravity can
be found in the GR literature. Refs. 2016JCAP…07..053A, 2017JCAP…11..004A and
2017arXiv170904703L present examples of implicit and explicit solvers that can
be used for going beyond the static limit in GR. Furthermore, as discussed in
Sec. 3, the equations associated with MG are not unique to gravity, but are
used also in different contexts such as geophysics, engineering (especially in
hydrodynamical applications) or medicine. Communication with these communities
is essential for making progress in this field.

##  5 Baryonic Simulations

The only difference that exists between simulating dark matter and simulating
baryons is that in the baryonic case, it is not possible to neglect
nongravitational interactions between the particles. This gives rise to a
plethora of thermodinamic and quantum processes between different elements of
the periodic table, transitions between them, transfer of energy between these
elements and photons, etc. The complexity of the solutions when including all
these effects increases disproportionally: A simulation that includes these
interactions and which has infinite resolution may end up producing humans
that may run simulations which may generate other humans that may run other
simulations and so on. Naturally, we do not have the computational resources
that are required to run such a simulation, so we need to make
simplifications. These simplifications are based on the idea of dividing the
solutions of the equations into different sets of scales. For the large
scales, the equations are solved self-consistently to follow the evolution of
dark matter and baryonic perturbations. The effects that occur at scales that
are smaller than the resolution of the simulation are treated analytically
following different processes with approximate models. As we are not
interested in studying the formation of humans in a cosmological context but
only galaxies, these sub-grid recipes have also a cut off in resolution.

The complexity of baryonic simulations increases as the field progresses. The
first baryonic simulations included only adiabatic gas and no sub-grid
physics. Later on, cooling of the gas was included, followed by recipes for
stellar formation, stellar evolution and different kinds of feedback (from
supernovae, black holes, etc). The state-of-the-art today includes magnetic
fields. Information about these kind of simulations within GR can be found in
the presentation of some of the latest simulations and references there in[78,
79, 80, 81]. The following sections describe how to include MG in these
simulations.

###  5.1 Including MG in hydrodynamical solvers

In order to understand the evolution of baryons in the nonlinear regime, we
need to reformulate the definition of cosmological simulation provided in Sec.
2 (Definition 2.2). In particular, we need to relax the assumption of
collisionless matter that is encoded in the energy-momentum tensor. Doing this
will prevent us from describing the fluid by following geodesics (Eq. (2)). So
to derive alternatives to these equations, we need to repeat the analysis that
was presented in Secs. 2.1 and 2.2, which resulted in Definition 2.3. I
provide a schematic of how this can be done below. More details on these
derivations can be found in the works that describe specific
implementations[82, 3, 83, 84].

Let us assume, for instance, a gravitational model described by the Lagrangian
(13). Minimizing the action with respect to the metric will give rise to the
following equation:

$G_{ab}=\frac{\delta\tilde{g}_{ab}}{\delta
g_{ab}}\tilde{T}_{ab}+T^{\phi}_{ab},$ (70)

where the two energy–momentum tensors are the usual Jordan frame matter tensor
and the Einstein frame energy–momentum tensor of the scalar field:

$\displaystyle\tilde{T}_{ab}$
$\displaystyle\propto\frac{1}{\sqrt{-g}}\frac{\delta
S_{M}}{\delta\tilde{g}_{ab}},$ (71) $\displaystyle T^{\phi}_{ab}$
$\displaystyle\propto\frac{1}{\sqrt{-g}}\frac{\delta S_{KG}}{\delta g_{ab}}.$
(72)

Taking the divergence of Eq. (70) gives a generalized equation of conservation
of energy. After projecting the result in the time and space directions, we
will end up with the usual continuity and Euler equations that can be found in
hydrodynamical codes plus MG terms that involve derivatives of
$T^{\phi}_{ab}$. In the Einstein frame, these new terms take the form of an
additional energy in the continuity equation and a fifth force in the Euler
equation and depend on the scalar field itself and its derivatives in space
and time. Thus, they can be obtained by solving the Klein–Gordon equation for
the field with the methods described in Sec. 3 or 4.

According to this analysis, implementing MG hydrodynamical solvers may be
straightforward: a MG hydrodynamical solver can be obtained simply by plugging
the solution of the Klein–Gordon equation that already exists in collisionless
MG codes into a standard hydrodynamical solver. However, additional complexity
may arise, which could lead to the need for more complex implementations. For
instance, in the case where the free function $f$ in Eq. (13) corresponds to a
standard kinetic term, the energy–momentum associated to the scalar field is

$T^{\phi}_{ab}=\nabla_{a}\phi\nabla_{b}\phi-\frac{1}{2}g^{cd}\nabla_{c}\phi\nabla_{d}\phi
g_{ab}-V(\phi)g_{ab},$ (73)

which is diagonal for a diagonal metric (i.e. when neglecting vector and
tensor modes). However, in the more general case in which $f$ includes cross
derivatives, this assumption may be broken, which may force us to take into
account nondiagonal components of the equations. Also the Klein–Gordon
equation for these models is sourced by the trace of the energy–momentum
tensor. In the case of baryons, there will be a pressure term sourcing the
equation, which is not included in existing simulations.

###  5.2 Including MG in sub-grid models and semi-analytics

The sub-grid models contain information about baryonic processes that occur at
the very small scales that lie below the resolution of the self-consistent
part of the simulation. These processes are certainly affected by the presence
of gravity and thus, depend on the assumed gravitational theory. For instance,
the formation of stars depends on how gas clouds fragment and collapse under
their own gravity. Different gravitational theories will give rise to
different evolution of the clouds. This in turn should give rise to a
different initial mass function, which is one of the ingredients in these sub-
grid models. In present simulations, these extremely complex gravitational
effects are included implicitly in the models (i.e. in most of the cases, the
gravitational constant does not appear explicitly in the equations). These
highly nonlinear phenomena are described by simple equations and condense
their complexity into a set of few free parameters. In the GR case, these
parameters are chosen by trial and error: Large number of test runs is
performed varying model parameters until a few observables such as the stellar
formation rate density, the luminosity function of galaxies, etc. agree with
the observed quantities. This process of tuning the parameters of the
simulations is very costly and done by large collaborations.

The generalization of the sub-grid models to MG is a two step process. First
it is necessary to re-derive the analytic expressions that describe different
aspects of the models. Second, the model parameters have to be recalculated.
In practice, constructing these models and tuning the simulations is an
extremely complex task, which can not be done by the small MG community. The
state-of-the-art simulations with MG do not include specifics of MG in the
definition of the models and instead use sub-grid models and parameters given
by the GR simulations that are usually run for comparison. Thus, existing
simulations are not yet self-consistent.

For models that include screening mechanisms, the approximation that is made
when not redefining the sub-grid models might be justified by the fact that
the impact of sub-grid physics is dominant in high density regions, where the
additional degrees of freedom should be screened. Examples of these kind of
simulations can be found in Refs. 2014MNRAS.442..176P and 2016PhRvD..93l3512H.
The impact of MG in some of the small scale problems that the $\Lambda$CDM
model historically faces, and whose proposed solutions within GR are usually
related to the impact of baryons, was presented in Refs.
corbett_moran_simulations, 2015MNRAS.453.1371M and 2016MNRAS.461.2490P.

A different approach usually applied in a GR context to the determination of
statistical properties of galaxies consists in taking into account baryonic
effects only in the post-processing phase of the simulations. In this case,
the simulations follow only the dark matter component. A set of recipes is
applied a posteriori on the resulting merger trees of dark matter halos to
determine the properties of the galaxies that form inside them. A change in
the underlying gravitational model, will affect the outcome of these
calculations through differences in both the merger trees and details on the
semi-analytic model itself. As in the case of sub-grid physics used in full
hydro simulations, applying GR semi-analytic codes to MG merger trees may not
give realistic results, specially in the MG models that include coupling.
Several quantities have to be redefined. A simple example is given by the
hydrostatic equilibrium, which is obtained by opposing pressure forces to
gravity forces and which is usually part of the semi-analytic models. In the
MG case, the fifth force should be included in the analysis. Another example
is the distribution of spin parameter of the dark matter halos. We know from
simulations that the distribution in MG is log-normal as in the GR case, but
the parameters (which are not free, but obtained from simulations) are
different. So a correct implementation of a semi-analytic model in MG will
require changing these and other similar quantities in the codes. Furthermore,
retuning of the free parameters of the models will be required. See Refs.
2013MNRAS.436.2672F and 2015MNRAS.452..978F for examples of MG gravity
implementations.

The question of whether the differences found between GR and MG galaxy
formation simulations are real or a consequence of poorly tuned baryonic
models is still open. For instance, differences were found in rotation curves
of galaxies in Ref. 2014MNRAS.442..176P. Are these differences based on
reality or simply on the fact that the simulations need re-tuning? Can we
recover the GR rotation curves in a MG context by retuning parameters? Only by
re-running the simulations we will be able to give answers to these questions.
This is extremely challenging, but by doing it we may open the door for
establishing new tests which may give us the definitive answer about gravity.
The risk is high, but certainly worth taking.

## Acknowledgements

CLL acknowledges support from STFC consolidated grant ST/L00075X/1 &
ST/P000541/1 and ERC grant ERC-StG-716532-PUNCA. CLL thanks Shaun Cole for
carefully reading the manuscript and anonymous referees and editor whose
comments greatly improve the quality of the paper.

## References

  * [1] J. Adamek, D. Daverio, R. Durrer and M. Kunz, J. Cosmology Astropart. Phys 7 (July 2016) 053, arXiv:1604.06065.
  * [2] J. Adamek, R. Durrer and M. Kunz, J. Cosmology Astropart. Phys 11 (November 2017) 004, arXiv:1707.06938.
  * [3] E. Puchwein, M. Baldi and V. Springel, MNRAS 436 (November 2013) 348, arXiv:1305.2418.
  * [4] C. Llinares, D. F. Mota and H. A. Winther, A&A 562 (February 2014) A78, arXiv:1307.6748.
  * [5] H. A. Winther, F. Schmidt, A. Barreira, C. Arnold, S. Bose, C. Llinares, M. Baldi, B. Falck, W. A. Hellwing, K. Koyama, B. Li, D. F. Mota, E. Puchwein, R. E. Smith and G.-B. Zhao, MNRAS 454 (December 2015) 4208, arXiv:1506.06384.
  * [6] A. Klypin, A. V. Macciò, R. Mainini and S. A. Bonometto, ApJ 599 (December 2003) 31, astro-ph/0303304.
  * [7] M. Baldi, V. Pettorino, G. Robbers and V. Springel, MNRAS 403 (April 2010) 1684, arXiv:0812.3901.
  * [8] A. V. Macciò, C. Quercellini, R. Mainini, L. Amendola and S. A. Bonometto, Phys. Rev. D 69 (June 2004) 123516, astro-ph/0309671.
  * [9] A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh and S. Pascoli, J. Cosmology Astropart. Phys 9 (September 2014) 031, arXiv:1408.1084.
  * [10] E. Carlesi, A. Knebe, G. Yepes, S. Gottlöber, J. B. Jiménez and A. L. Maroto, MNRAS 424 (July 2012) 699, arXiv:1205.1695.
  * [11] H. Oyaizu, Phys. Rev. D 78 (December 2008) 123523, arXiv:0807.2449.
  * [12] B. Li, G.-B. Zhao, R. Teyssier and K. Koyama, J. Cosmology Astropart. Phys 1 (January 2012) 051, arXiv:1110.1379 [astro-ph.CO].
  * [13] P. Brax, A.-C. Davis, B. Li, H. A. Winther and G.-B. Zhao, J. Cosmology Astropart. Phys 10 (October 2012) 002, arXiv:1206.3568.
  * [14] A.-C. Davis, B. Li, D. F. Mota and H. A. Winther, ApJ 748 (March 2012) 61, arXiv:1108.3081 [astro-ph.CO].
  * [15] F. Schmidt, Phys. Rev. D 80 (August 2009) 043001, arXiv:0905.0858 [astro-ph.CO].
  * [16] J. Khoury and M. Wyman, Phys. Rev. D 80 (September 2009) 064023, arXiv:0903.1292 [astro-ph.CO].
  * [17] I. Laszlo and R. Bean, Phys. Rev. D 77 (January 2008) 024048, arXiv:0709.0307.
  * [18] K. C. Chan and R. Scoccimarro, Phys. Rev. D 80 (November 2009) 104005, arXiv:0906.4548 [astro-ph.CO].
  * [19] B. Li, G.-B. Zhao and K. Koyama, J. Cosmology Astropart. Phys 5 (May 2013) 023, arXiv:1303.0008 [astro-ph.CO].
  * [20] A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh and S. Pascoli, J. Cosmology Astropart. Phys 10 (October 2013) 027, arXiv:1306.3219.
  * [21] B. Li, A. Barreira, C. M. Baugh, W. A. Hellwing, K. Koyama, S. Pascoli and G.-B. Zhao, J. Cosmology Astropart. Phys 11 (November 2013) 012, arXiv:1308.3491 [astro-ph.CO].
  * [22] C. Llinares and D. F. Mota, Physical Review Letters 110 (April 2013) 161101, arXiv:1302.1774 [astro-ph.CO].
  * [23] C. Llinares and D. F. Mota, Phys. Rev. D 89 (April 2014) 084023, arXiv:1312.6016 [astro-ph.CO].
  * [24] S. Bose, W. A. Hellwing and B. Li, J. Cosmology Astropart. Phys 2 (February 2015) 034, arXiv:1411.6128.
  * [25] R. Hagala, C. Llinares and D. F. Mota, A&A 585 (January 2016) A37, arXiv:1504.07142.
  * [26] A. Nusser, MNRAS 331 (April 2002) 909, astro-ph/0109016.
  * [27] A. Knebe and B. K. Gibson, MNRAS 347 (February 2004) 1055, astro-ph/0303222.
  * [28] R. Brada and M. Milgrom, ApJ 519 (July 1999) 590, astro-ph/9811013.
  * [29] O. Tiret and F. Combes, A&A 464 (March 2007) 517, astro-ph/0701011.
  * [30] C. Llinares, A. Knebe and H. Zhao, MNRAS 391 (December 2008) 1778, arXiv:0809.2899.
  * [31] C. Llinares, On the linear and non-linear evolution of dust density perturbations with MOND, PhD thesis, Rijksuniversiteit Groningen (2011).
  * [32] P. Londrillo and C. Nipoti, Memorie della Societa Astronomica Italiana Supplementi 13 (2009) 89, arXiv:0803.4456.
  * [33] P. Londrillo and C. Nipoti, N-MODY: A Code for Collisionless N-body Simulations in Modified Newtonian Dynamics Astrophysics Source Code Library (February, 2011).
  * [34] G. W. Angus, K. J. van der Heyden, B. Famaey, G. Gentile, S. S. McGaugh and W. J. G. de Blok, MNRAS 421 (April 2012) 2598, arXiv:1201.3185.
  * [35] G. N. Candlish, R. Smith and M. Fellhauer, MNRAS 446 (January 2015) 1060, arXiv:1410.3844.
  * [36] F. Lüghausen, B. Famaey and P. Kroupa, Canadian Journal of Physics 93 (February 2015) 232, arXiv:1405.5963.
  * [37] C. Llinares, ArXiv e-prints (September 2017) arXiv:1709.04703.
  * [38] R. W. Hockney and J. W. Eastwood, Computer simulation using particles 1988\.
  * [39] S. M. Carroll, Spacetime and geometry. An introduction to general relativity 2004\.
  * [40] H. Martel and P. R. Shapiro, MNRAS 297 (June 1998) 467, astro-ph/9710119.
  * [41] N. E. Chisari and M. Zaldarriaga, Phys. Rev. D 83 (June 2011) 123505, arXiv:1101.3555 [astro-ph.CO].
  * [42] S. R. Green and R. M. Wald, Phys. Rev. D 85 (March 2012) 063512, arXiv:1111.2997 [gr-qc].
  * [43] C. Fidler, C. Rampf, T. Tram, R. Crittenden, K. Koyama and D. Wands, Phys. Rev. D 92 (December 2015) 123517, arXiv:1505.04756.
  * [44] C. Fidler, T. Tram, C. Rampf, R. Crittenden, K. Koyama and D. Wands, J. Cosmology Astropart. Phys 12 (December 2017) 022, arXiv:1708.07769.
  * [45] E. Bertschinger, ARA&A 36 (1998) 599.
  * [46] J. S. Bagla, Current Science 88 (April 2005) 1088, astro-ph/0411043.
  * [47] K. Dolag, S. Borgani, S. Schindler, A. Diaferio and A. M. Bykov, Space Sci. Rev. 134 (February 2008) 229, arXiv:0801.1023.
  * [48] S. Borgani and A. Kravtsov, Advanced Science Letters 4 (February 2011) 204, arXiv:0906.4370 [astro-ph.CO].
  * [49] M. Zumalacárregui, T. S. Koivisto and D. F. Mota, Phys. Rev. D 87 (April 2013) 083010, arXiv:1210.8016 [astro-ph.CO].
  * [50] R. P. Woodard, Foundations of Physics 44 (February 2014) 213, arXiv:1401.0254.
  * [51] M. Maggiore and M. Mancarella, Phys. Rev. D 90 (July 2014) 023005, arXiv:1402.0448 [hep-th].
  * [52] Y. Dirian, S. Foffa, N. Khosravi, M. Kunz and M. Maggiore, J. Cosmology Astropart. Phys 6 (June 2014) 033, arXiv:1403.6068.
  * [53] W. Hu and I. Sawicki, Phys. Rev. D 76 (September 2007) 064004, arXiv:0705.1158.
  * [54] P. Brax, C. van de Bruck, A.-C. Davis and D. J. Shaw, Phys. Rev. D 78 (November 2008) 104021, arXiv:0806.3415.
  * [55] S. Bose, B. Li, A. Barreira, J.-h. He, W. A. Hellwing, K. Koyama, C. Llinares and G.-B. Zhao, J. Cosmology Astropart. Phys 2 (February 2017) 050, arXiv:1611.09375.
  * [56] G. Dvali, G. Gabadadze and M. Porrati, Physics Letters B 485 (July 2000) 208, hep-th/0005016.
  * [57] F. Schmidt, Phys. Rev. D 80 (December 2009) 123003, arXiv:0910.0235 [astro-ph.CO].
  * [58] M. Milgrom, ApJ 270 (July 1983) 365.
  * [59] J. Bekenstein and M. Milgrom, ApJ 286 (November 1984) 7.
  * [60] M. Milgrom, MNRAS 403 (April 2010) 886, arXiv:0911.5464.
  * [61] M. Milgrom, Phys. Rev. D 80 (December 2009) 123536, arXiv:0912.0790 [gr-qc].
  * [62] A. Brandt, Math. of Comp. 31 (April 1977) 333.
  * [63] P. Wesseling, An Introduction to Multigrid Methods (John Wiley and Sons Inc, |c1992, 2nd ed., 1992).
  * [64] U. Trottenberg, C. Oosterlee and A. Scholler, Multigrid (Academic Press, 2000).
  * [65] A. Brandt and O. Livne, Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics, Revised EditionClassics in Applied Mathematics, Classics in Applied Mathematics (Society for Industrial and Applied Mathematics, 2011).
  * [66] A. Barreira, B. Li, C. M. Baugh and S. Pascoli, J. Cosmology Astropart. Phys 11 (November 2013) 056, arXiv:1308.3699.
  * [67] H. A. Winther and P. G. Ferreira, Phys. Rev. D 92 (September 2015) 064005, arXiv:1505.03539 [gr-qc].
  * [68] A. Barreira, S. Bose and B. Li, J. Cosmology Astropart. Phys 12 (December 2015) 059, arXiv:1511.08200.
  * [69] H. Weinberger, A First Course in Partial Differential Equations with Complex Variables and Transform MethodsBlaisdell book in pure and applied mathematics, Blaisdell book in pure and applied mathematics (Dover Publications, 1995).
  * [70] R. Agarwal and D. O’Regan, Ordinary and Partial Differential EquationsUniversitext, Universitext (Springer-Verlag New York, 2009).
  * [71] J. Noller, F. von Braun-Bates and P. G. Ferreira, Phys. Rev. D 89 (January 2014) 023521, arXiv:1310.3266 [astro-ph.CO].
  * [72] L. M. Widrow, Phys. Rev. D 40 (August 1989) 1002.
  * [73] W. H. Press, B. S. Ryden and D. N. Spergel, ApJ 347 (December 1989) 590.
  * [74] R. Hagala, C. Llinares and D. F. Mota, Physical Review Letters 118 (March 2017) 101301, arXiv:1607.02600.
  * [75] C. Llinares and L. Pogosian, Phys. Rev. D 90 (December 2014) 124041, arXiv:1410.2857.
  * [76] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics - A practical introduction - 2nd edition (Springer, Berlin, 1999).
  * [77] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3 edn. (Cambridge University Press, New York, NY, USA, 2007).
  * [78] M. Vogelsberger, S. Genel, D. Sijacki, P. Torrey, V. Springel and L. Hernquist, MNRAS 436 (December 2013) 3031, arXiv:1305.2913.
  * [79] M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu, G. Snyder, D. Nelson and L. Hernquist, MNRAS 444 (October 2014) 1518, arXiv:1405.2921.
  * [80] J. Schaye, R. A. Crain, R. G. Bower, M. Furlong, M. Schaller, T. Theuns, C. Dalla Vecchia, C. S. Frenk, I. G. McCarthy, J. C. Helly, A. Jenkins, Y. M. Rosas-Guevara, S. D. M. White, M. Baes, C. M. Booth, P. Camps, J. F. Navarro, Y. Qu, A. Rahmati, T. Sawala, P. A. Thomas and J. Trayford, MNRAS 446 (January 2015) 521, arXiv:1407.7040.
  * [81] A. Pillepich, V. Springel, D. Nelson, S. Genel, J. Naiman, R. Pakmor, L. Hernquist, P. Torrey, M. Vogelsberger, R. Weinberger and F. Marinacci, MNRAS 473 (January 2018) 4077, arXiv:1703.02970.
  * [82] L. Casarini, A. V. Macciò, S. A. Bonometto and G. S. Stinson, MNRAS 412 (April 2011) 911, arXiv:1005.4683.
  * [83] A. Hammami, C. Llinares, D. F. Mota and H. A. Winther, MNRAS 449 (June 2015) 3635, arXiv:1503.02004.
  * [84] J.-h. He and B. Li, Phys. Rev. D 93 (June 2016) 123512, arXiv:1508.07350.
  * [85] C. Penzo, A. V. Macciò, L. Casarini, G. S. Stinson and J. Wadsley, MNRAS 442 (July 2014) 176, arXiv:1401.3338.
  * [86] C. C. Moran, R. Teyssier and B. Li, Monthly Notices of the Royal Astronomical Society 448 (2015) 307.
  * [87] A. V. Macciò, R. Mainini, C. Penzo and S. A. Bonometto, MNRAS 453 (October 2015) 1371, arXiv:1503.07867.
  * [88] C. Penzo, A. V. Macciò, M. Baldi, L. Casarini, J. Oñorbe and A. A. Dutton, MNRAS 461 (September 2016) 2490, arXiv:1504.07243.
  * [89] F. Fontanot, E. Puchwein, V. Springel and D. Bianchi, MNRAS 436 (December 2013) 2672, arXiv:1307.5065.
  * [90] F. Fontanot, M. Baldi, V. Springel and D. Bianchi, MNRAS 452 (September 2015) 978, arXiv:1505.02770.