Ecosmog: An Efficient Code for Simulating Modified Gravity
Abstract
We introduce a new code, ECOSMOG, to run body simulations for a wide class of modified gravity and dynamical dark energy theories. These theories generally have one or more new dynamical degrees of freedom, the dynamics of which are governed by their (usually rather nonlinear) equations of motion. Solving these nonlinear equations has been a great challenge in cosmology. Our code is based on the RAMSES code, which solves the Poisson equation on adaptively refined meshes to gain high resolutions in the highdensity regions. We have added a solver for the extra degree(s) of freedom and performed numerous tests for the gravity model as an example to show its reliability. We find that much higher efficiency could be achieved compared with other existing mesh/gridbased codes thanks to two new features of the present code: (1) the efficient parallelisation and (2) the usage of the multigrid relaxation to solve the extra equation(s) on both the regular domain grid and refinements, giving much faster convergence even under much more stringent convergence criteria. This code is designed for performing highaccuracy, highresolution and largevolume cosmological simulations for modified gravity and general dark energy theories, which can be utilised to test gravity and the dark energy hypothesis using the upcoming and future deep and highresolution galaxy surveys.
pacs:
I Introduction
The accelerating expansion of our Universe is one of the most challenging questions in modern physics cst2006 (). After more than a decade of attempts in constructing consistent and natural theories for it, we are still nowhere near a definite conclusion. Broadly speaking, the theoretical models developed so far roughly fall into two categories: those involving some exotic matter species, the socalled dark energy, which usually have nontrivial dynamics, and those involving modifications to the Einstein gravity on certain (usually large) scales. Examples of the dynamical dark energy include the quintessence zws1999 (), kessence ams2000, coupled quintessence a2000 (), chameleon field kw2004 (); ms2007 (), symmetron field hk2010 () etc., and some examples of modified gravity models are the gravity cddett2005 (), scalartensor theory pb1999 () and DGP model dgp2000 () (for some recent review see, e.g., k2008 (); dm2008 (); sf2008 ()).
From a practical point of view (i.e., regardless of the naturalness or finetuning considerations), there are two important issues faced by any theoretical model: consistency and degeneracy.
In order to ensure the consistency of the mode, a given model should not violate any existing observational constraints. The new degrees of freedom, which are supposed to drive the accelerating expansion only on very large scales, quite often produce unwanted side effects. Consider the coupled quintessence as an example; if it couples to normal matter species (baryons etc.), it can mediate a fifth force that is strongly constrained by experiments. This problem can certainly be avoided by assuming that the quintessence field only couples to dark matter which we cannot directly observe, but more interesting theoretical mechanisms to avoid this problem have been designed in recent years. The leading example is the chameleon mechanism kw2004 (); ms2007 (); lb2007 (). By this mechanism, the fifth force is suppressed to the undetectable level in regions with high matter density, while in lowdensity environments such as galaxies, galaxy clusters and cosmological voids, the fifth force is unsuppressed and can be as strong as gravity. The rapid changes of the behaviour of the fifth force across different regions naturally imply that the equation(s) of motion governing the dynamics of the new degree(s) of freedom should be very nonlinear. This adds the complexity to the model but we cannot avoid this to ensure the consistency of the model.
Different theoretical models can behave very similarly to each other and to the standard CDM paradigm, which makes it difficult to distinguish amongst them using observational data. The (apparent) existence of the degeneracy usually indicates the incompleteness of our theoretical understanding of models. For instance, it is very easy to design a scalar field model which produces almost identical cosmic expansion history as CDM, but this is merely because in the background cosmology, any detailed structure of the Universe is disregarded; when we look at the linear perturbation evolutions, the apparent degeneracy is often broken. However, there are situations where the degeneracy cannot be broken even with linear perturbation analysis (such as the chameleon theory and gravity in which the length scales on which linear perturbation analyses apply are bigger than the range of, and so not affected by, the fifth force). In such situations, a full study of the nonlinear structure formation and evolution needs to be performed to break the degeneracy.
Therefore, it is essential to study nonlinea structure formation to ensure the consistency of the model and break the degeneracies between various models. Nonlinear structure formation is too complex to understand analytically and therefore requires numerical simulations. Numerical simulations of the cosmic structure evolution on small (such as galactic or cluster) scales can not only open a new arena of using consistency with observations to constrain models, but also hopefully break the degeneracy between different models (as our cosmological simulations show below). This is particularly exciting considering that highquality observational data will keep coming in the next decades to improve our theoretical understanding.
Numerical simulations are done by numerical codes. In most cases the extra dynamical degree(s) of freedom is more accurately treated as a field and its value is required on a number of space points. For this purpose we need the numerical code that can solve the field value on meshes covering (parts of) the simulation box. There are two such codes known to us to date: one is that of Oyaizu o2008 (), which has been applied to gravity olh2008 (); sloh2009 () and the DGP model f2009 (); the other is a modified MLAPM code mlapm () written by one of the authors lz2009 () and which has been applied to chameleon theories lz2010 (), gravity zlk2011 (), coupled scalar field theories lb2011 (); lmb2011a (), scalartensor theory lmb2011b (), dilaton model bbdls2011 (), symmetron model dlmw2011 (). Both have shortcomings: the Oyaizu’s code does not support adaptive refinements (thus easier to implement) and highresolution simulations are not practical; while the MLAPM code does support adaptive refinements, it is not parallelised and not practical for simulations with very big volumes and high mass resolutions. Furthermore, the equations on MLAPM refinements are solved on a onelevel grid, which is rather inefficient.
The aim of this paper is to introduce a new code ECOSMOG, which overcomes the shortcomings of the previous codes. The new code is based on RAMSES ramses (). It is efficiently parallelised, supports adaptive refinements and solves the equations using multigrid method on the refinements (for a more detailed comparison of the three codes the readers can refer to the table I). As a working example, we use the new code to run a number of test simulations for the gravity, which is one of the most challenging models to simulate because of the high nonlinearity of its equations. As will be shown below, the code works very well for the model, and we certainly expect it to be straightforward to implement equations in other models to our code.
The organisation of this paper is as follows. In § II we briefly introduce the gravity model. In § III we introduce the supercomoving code unit using a different form and list the body Poisson and equations to simplify the numerics. § IV then makes discrete versions of these equations to be implemented in our code. § V, we show details on how the numerical implementation is performed and discuss several important issues, which is the core section of this paper. Next, a long section, § VI, contains the results of eleven tests of the code. These tests check the correctness, efficiency and consistency of different aspects of the code and give us confidence about its reliability. Finally we compare the present code with other meshbased codes, summarise and conclude in § VII.
This is a paper to explain the code and physical interpretations of the results will be presented in future publications.
Ii A Test Case: The Gravity
One can alter the Einstein gravity in such a way that it gives rise the cosmic acceleration without introducing dark energy. One example along this line is the socalled gravity, in which the Ricci scalar in the EinsteinHilbert action is generalised to a function of (see e.g., sf2008 () for a review and references therein). In gravity, the structure formation is governed by the following two equations,
(1)  
(2) 
where denotes the gravitational potential, is the extra scalar degree of freedom, dubbed scalaron, , and the quantities with overbar take the background values. The symbol is the three dimensional gradient operator, and is the scale factor. These two coupled Poissonlike equations are much difficult to solve than the single Poisson equation in General Relativity (GR), , which is linear. From Eqs (1) and (2), we can see that gravity in can be enhanced depending on the local environment – in underdense regions, the term in Eqs (1) vanishes thus the two equations decouple, making gravity simply enhanced by a factor of . However in the dense region, in Eq (2) is negligible, yielding , which means that GR is locally restored. This is the chameleon mechanism applied in the gravity models, which is important for the cosmological viability of the latter. The presence of the chameleon effect indicates that Eqs (1) and (2) are highly nonlinear, making it challenging to numerically solve them in the simulation process.
An gravity model is fully specified by the functional form of , and here we shall adopt the model proposed by Hu & Sawicki hs2007 (), which takes the form
(3) 
where are model parameters, and
(4) 
with being the present fractional matter density and the current Hubble expansion rate.
It can then be shown that
(5) 
Because
(6) 
where overbars are used for the background quantities, to match a CDM background expansion we have . With and , we find , which means that can be approximated as
(7) 
Because is the actual quantity that enters the body equations (see below), we find that only the two parameters and , are needed in our simulations. Another (independent) parameter which is the present background value of , can be obtained from and is often used to give people some idea about the size of .
Iii The body Equations
In this section we shall introduce the convention of our code units, and list the body Poisson and equations, the latter being written in the socalled quasistatic limit so that terms involving time derivatives will be dropped. The body equations can all be found elsewhere, though perhaps of slightly different forms; consequently, this section shall be kept short and only serves for completeness.
iii.1 The Code Units
The code unit used in the RAMSES code and its modification developed here is based on (but not exactly) the supercomoving coordinates of ms1998 (). It can be summarised as follows (where tilded quantities are expressed in the code unit):
In the above is the comoving coordinate, is the critical density today, the fractional energy density for matter today, the particle velocity, the gravitational potential and the speed of light. In addition, is the size of the simulation box in unit of Mpc and the Hubble expansion rate today in unit of km/s/Mpc. Note that with these conventions the average matter density is . All the newlydefined quantities are dimensionless.
In the gravity equation we also have a new degree of freedom and in the code unit we will use instead. As we shall see below, this is to make sure that has an equivalent status as : the latter is the Newtonian potential and determines the total (modified) gravitational force, while the former is related to the potential of the extra force and determines the modification to the standard Einsteinian gravity.
iii.2 The Modified Poisson Equation
From above we see that the modified Poisson equation is
(8) 
where and
(9)  
(10) 
Using the code units defined above, the equation becomes
(11)  
iii.3 The Equation
The equation of motion for ,
(12) 
becomes in the code units,
(13)  
Comparing this equation with that for the modified Poisson equation, we find that and has the same code unit, and that is why we have used instead of .
As in zlk2011 (), we will not solve directly as it may change too quickly from one space point to another and can cause divergence problems in the numerical solutions. Instead, we will solve a new variable , defined through . will be of order unity across the space, and so has better convergence properties.
Iv The Discretised Equations
From here on we will only use variables expressed using the code unit, and the tildes will be dropped for simplicity.
Clearly, to put the above equations into the body code one must discretise them. For the Poisson equation, it is linear as along as the source term (the righthand side) is provided, so the discretisation is fairly straightforward:
(14)  
where for example is the value of in the grid cell with index . The discrete version of the nonlinear equation of motion, in contrast, is more tedious zlk2011 ():
(15) 
Notice that in the above equations we have used to avoid confusion in notations and to make it clear that we are using the density contrast. and
V The body Algorithms
The detailed implementation of the body solvers in the RAMSES code is indeed quite different from that in the MLAPM code lz2009 (); lb2011 (). As a result, certain corrections need to be made to the above discrete equations before implementing them.
Both codes adaptively refine the grid in highdensity regions and solve the Poisson and equations on the refinement to get better spatial resolutions. The refining is performed on a gridbygrid basis so that the final refinements typically have irregular shapes roughly matching the isodensity surfaces. Particles in regions underlying the refinements are linked to the latter, where their positions and momenta are updated.
An important difference between the two codes is how the equations are solved on the refinements. MLAPM solves the equations on single level of the considered refinement only, performing GaussSeidel relaxations till the convergence of the solution is obtained. In RAMSES, however, for each adaptive mesh refinement (AMR) level there is a series of separate, coarser multigrid levels on which the GaussSeidel relaxation is used to accelerate the convergence.
The treatments of boundary conditions on refinements are also different. In MLAPM, the outermost cells of a refinement are taken as its physical boundary, thus the field values wherein are set using interpolation from the coarserlevels solutions. In RAMSES, the use of multilevel requires the boundary conditions to be set consistently on the different levels, and this is achieved using an elegant mask scheme: on the AMR level, which is the finest level in the multigrid series, the active cells, inside which the field values need to be updated during the relaxation process, are given a mask value of while the other cells are assigned a mask value . The boundary of the computational domain is defined to be where the linearlyinterpolated mask value vanishes, or equally the boundary of the active AMR cells. The value of the field on this physical boundary is computed at the beginning of the relaxation process and kept unchanged after that until a converged solution in the active cells are obtained. The mask values in the coarser multigrid cells are restricted from their finerlevel values, while the physical boundary on this level is still defined as where the mask hits zero, which ensures the boundaries on different levels are consistent (for more details and tests in the classical GR case see gt2011 ()).
Moreover, the treatments of the boundary conditions for linear and nonlinear elliptical PDEs are also different in RAMSES, as we shall see below.
v.1 The Modified Poisson Equation
v.1.1 The BCmodified Equation
Unlike the MLAPM code, where the refinement boundary consists of the outermost cells on that refinement and so the boundary cells are always there, in RAMSES we might well face the situation that an outermost active cell has no neighbouring cells in some directions on the save level (possibly because a neighbouring coarse cell has not been refined), while we do need the field values in those cells to interpolate and compute the boundary conditions. Such cells are called ghost cells: they do not exist in the code data structure but we do need them.
For the linear elliptical PDEs, such as the (modified) Poisson equation, there is a simple way to avoid storing information about the ghost cells. Consider the discrete Poisson equation and now suppose that the cell is a ghost, which has mask value ^{1}^{1}1If this is a ghost cell then its mask value will be but here we try to be general in our description.. We realise that the physical boundary, namely where the mask value crosses zero, will be a distance
from the (ghost!) cell and
from the cell, where is the cell length, or equally the distance between the centres of these two cells. Using linear interpolation, the boundary value of will then be
which gives
(16) 
Note that though the cell is a ghost cell, the value of can be computed, at the beginning of the relaxation iterations, from its father (coarser) cell which does exist. Then is computed and fixed (because it is the boundary value!). During the subsequent relaxation iterations, the value of changes and so does that of but not . As we don’t have the cell, there is nowhere to store the updated values of , and so we choose to not use it at all, replacing in the discrete equation by and using Eq. (16). We then get a boundarycondition (BC)modified equation
(17)  
Note that because is fixed, we could simply move it to the righthand side of the equation. The BCmodified righthand side is then only needed to be computed once, before entering the relaxation iterations. There is no need to store either or . Note that the lefthand side is also modified and one must be careful here.
v.1.2 The Multigrid Implementation
Now consider the multigrid algorithm. For simplicity let us rewrite the above BCmodified equation as
(18) 
in which the superscript means we are on the level with cell size . Note that is a linear operator, which is important here.
Suppose after a number of GaussSeidel iterations (presmoothing), we get an approximate solution , then
(19) 
in the active cells, and the difference
(20) 
is called the residual. We could define and rewrite the equation as
(21) 
and coarsify this equation as
(22) 
and solve it on the coarser level to accelerate the convergence. Here, the superscript indicates that we are on the coarser level on which the cell size is , and is the restriction operator.
To solve this coarsified equation, we need the boundary conditions for on the coarser level and, as in Eq. (17), the coarserlevel equation will be modified to include correction terms involving , which is the boundary on the coarser level. However, on the boundary vanishes identically and it turns out that the coarserlevel version of Eq. (17) is simpler:
(23) 
assuming that the coarserlevel cell is a ghost cell. Note that the righthand side of the coarsified equation is not BC modified since , but the lefthand side is modified by the lack of neighbouring cells.
Here we have only described the algorithm for two levels, but the generalisation to more levels is trivial.
Finally some comments should be made on the restriction operation : if a given fine cell is masked^{2}^{2}2Here ”masked” means having a nonpositive mask value. then it has no contribution to the coarserlevel residual, because the field value in this cell is unmodified by the relaxation and so considered to be exact. At the same time, we shall not restrict residuals into coarse cells which are masked, because this is unnecessary.
v.1.3 The Prolongation
Once an approximate coarserlevel solution to is obtained, say , one can correct the finelevel solution using
(24) 
where is the prolongation operator, to find (expected) more accurate solutions on the fine level.
In the prolongation process, a fine cell receives contribution from its eight neighbouring coarser cells. Of these eight coarse cells:

one contains the fine cell and is given a weight of ,

three have one common face with the fine cell and are given a weight of each,

three have one common edge with the fine cell and are given a weight of each,

one has a common vertex with the fine cell and is given a weight of .
Of course, these weights sum to unity.
Note that if a given fine cell is masked, then it is outside the active computational domain and corrections are not necessary for it. If the coarse cell is not a valid multigrid cell (it could be a valid AMR cell, though), then its contribution to the finecell correction is zero.
v.2 The Equation
Having recalled the multigrid algorithm for the linear (modified) Poisson equation, first presented in gt2011 (), we can now have a look at the nonlinear equation. The nonlinearity introduces certain complications compared to the linear case, which should be taken care of in the numerical implementations.
v.2.1 The BCModified Equation
As in the linear case, let us suppose that the cell on the finest level is a ghost cell. Then
(25) 
with the boundary value of . Note that this necesarilly means that the equality
(26) 
does not hold, because is nonlinear in – this point is very important in order to make consistent BCmodified equation. Instead, we shall use with given by Eq. (25).
As before, is computed and then used to obtain before entering the relaxation iterations. In the subsequent iterations, is updated and so is , but not .
Removing the and in Eq. (IV) using the above two relations, we arrive at the BC modified equation:
(27)  
Now we could see clearly where the additional complexity appears. Remember that in the linear case we don’t have to store because it only appears on the BCmodified righthand side of the equation. Here, on the other hand, also appears in the term such that its value is needed in each iteration during which is updated. This implies that, for each outermost active cell, we should store for its neighbouring ghost cell(s).
If we define the the lefthand side of Eq. (27) as for simplicity, then it can be easily shown that
(28)  
Note that when and , the above equations reduce to those for a regular grid with periodic boundary conditions as expected. The relaxation on this level is then done through
(29) 
v.2.2 The Multigrid Implementation
Let us write the BCmodified equation as
(30) 
on the fine level, where is the nonlinear operator acting on defined above. Note that here , but we shall still keep it for a reason which will become clear below.
After a number of presmoothing iterations, one finds an approximate solution on the fine level, which gives
(31) 
Consider the new equation
(32) 
After coarsifying and rearranging, we obtain the coarserlevel equation as
(33) 
After a number of relaxation iterations on the coarser level, we find an approximate solution , and the finelevel solution can be corrected as
(34) 
where is the prolongation operator.
Different from the linear case, here on the coarser level we are not solving an equation for the correction to the field (remember in that case we solved on the coarser level) but again itself. Therefore Eqs. (27, 28) could be applied to the coarse level as well, if we

change to , and

find correct coarserlevel boundary values for , .
The first point is fairly straightforward, while the second needs some further analysis. One must first find the physical boundary on the coarser level which, as explained above, is where the coarserlevel mask crosses zero. This is easy because is computed by restricting the finecell mask values. The values of in the boundary coarse cells are taken as those in the corresponding AMR cells.
Finally some comments should be made on the restriction operations and . will be treated similarly as in the linear case for the same reasons explained there.
For , even though a given fine cell is masked it still has contribution, because otherwise the computed for the coarser cell will be incorrect. But if a coarser cell is masked we will not compute for it, since this will be unused anyway, as in the case of .
The truncation error can be estimated as
(35) 
and similarly for other levels. This could provide a stopping (convergence) criterion for the multigrid iteration, which in our case is
(36) 
where is a predefined constant.
v.2.3 The Prolongation
As mentioned above, in the nonlinear case we have to prolongate the quantity , but not the residual itself. The prolongated result can then be used to correct the finelevel solutions.
The prolongation here is done in the same way as that for the residual of the linear (modified) Poisson equation, to be consistent. Details shall not be presented here.
Vi Code Tests
In this section we shall give the results of some tests we have performed to show that the above algorithm and implementation work correctly and efficiently.
vi.1 Equation Solver On Domain Grid
The most important part in the modified RAMSES code, which is the topic of the last two sections, is the equation of motion for gravity or the new degree of freedom(s) in other theories. Here we have performed a range of tests for it with different configurations of the matter density field.
vi.1.1 Homogeneous Matter Density Field
In a universe with homogeneous density, we know that the quantity is homogeneous and it exactly takes its background value , namely
(37) 
everywhere. Therefore, as the simplest test of the equation solver, one has to show that in such a homogeneous field, given some random guess of on the cells of the simulation mesh, after a reasonable number of GaussSeidel relaxation sweeps, the solution converges to the above background value. Such simple test have been used previously in bbdls2011 (); dlmw2011 () to show that the MLAPM solver for extra degrees of freedom works well.
We have performed this test for and a value of corresponding to . The result is shown in Fig. 1, where we plot the values of in the cells in direction, before and after the GaussSeidel relaxation, for two different random initial guesses. We can see that the final solution agrees with the analytical result (the horizontal line) very well (see figure caption for more details).
vi.1.2 Point Mass
As a second test of the RAMSES equation solver, let us consider the solution of around a point mass at the origin, for which case we have an analytical solution which is accurate except for the regions very close to the mass. Such a test has been used previously in o2008 (); bbdls2011 ().
Following o2008 (), we construct the pointmass density field as
(38) 
in which are respectively the cell indices in direction. In the test we have used a cubic box with size Mpc and 128 grid cells in each direction. The other physical parameters are chosen as and we have used three values of corresponding to .
Meanwhile, the analytical solution can be obtained approximately by solving the equation
(39) 
in which the effective mass of the scalar degree of freedom , , is given by
(40) 
Fig. 2 shows the comparison between the numerical solutions to along the axis (symbols) and analytical solutions (solid curves), and we can see that the two agree very well for all three models used in this test.
vi.1.3 Sine Density Field
As our third test, let us consider the sine density field introduced in o2008 (), which (after some modification to account for the code units) is given by
(41)  
in which is rescaled such that . Notice that we consider only dependence, which is equivalent to a onedimensional configuration. The solution to this density field can be analytically worked out to be,
(42) 
Fig. 3 shows the test results for the sine density field given above, with and three values of corresponding to respectively. It can be seen that the numerical solutions (symbols) agree with the analytical solutions (solid curves) very well.
vi.1.4 Gaussian Density Field
The last test of the equation solver on the domain grid uses a Gaussian type density configuration. Again, here we only consider onedimensional case, where the density field is specified by
(43)  
where again has been scaled to code units so that , , are numerical constants which respectively specify the width and height of the density field, which obviously peaks at , and is defined by
(44) 
Note that such a density field is not exactly periodic at the edges of the simulation box, but given that is small enough, at the box edges and periodic boundary conditions are approximately satisfied.
The solution to can then be obtained analytically and is
(45) 
which clearly shows that when could be made very small at while at or goes to its background value.
We have implemented Eq. (43) into our numerical code and the numerical solutions for are shown in Fig. 4. For simplicity we only show the results for , but for three different values of and (the open circles, open triangles and open squares in Fig. 4 respectively). We can see that the numerical results agree with the analytical solution Eq. (45) very well.
vi.2 Equation Solver On Refinements
The above four tests show that our solver for the equation actually works accurately on the regular domain grid, but this is not sufficient for the code test since the equation is also solved on refinements where, as we have seen, the equations take different forms due the complex treatment of the boundary conditions. It is therefore essential to test the equation solver on the refinements as well, as we shall do in this subsection.
vi.2.1 Twolevel Gaussian Density Field
The Gaussiantype density configuration could provide a good check of the multilevel equation solver because the density peak can be made arbitrarily high by adjusting the parameter . Near its peak, the density field changes rapidly and a higher spatial resolution is needed to compute accurately. Consider the case where the regular domain grid is refined only once, in regions where the density value exceeds some certain threshold (we shall call this a twolevel problem, and in the present numeric example the coarse and fine levels are respectively levels 8 and 9). The density values in both the coarse and refined cells are given by Eq. (43), while the values of at the finelevel boundaries are computed from interpolation of those in the nearby coarselevel cells, as discussed above.
Fig. 5 shows the numerical values of on both levels in the region covered by the refinement. We show the results for four different values of (, , from top to bottom), and for each the results from the coarse and fine levels are denoted respectively by filled squares and empty circles. For comparison we have also plotted the analytical results Eq. (45) as solid curves. As we can see, both finelevel and coarselevel results are virtually indistinguishable from the exact solution by eye.
This does not mean that the refinement is unnecessary however, because, as shown in Fig. 5, the fine level has more data points and could probe regions closer to the extreme value of , which corresponds to the highdensity region where high resolution is needed.