SelfConsistent Cosmological Simulations of DGP Braneworld Gravity
Abstract
We perform cosmological Nbody simulations of the DvaliGabadadzePorrati braneworld model, by solving the full nonlinear equations of motion for the scalar degree of freedom in this model, the brane bending mode. While coupling universally to matter, the branebending mode has selfinteractions that become important as soon as the density field becomes nonlinear. These selfinteractions lead to a suppression of the field in highdensity environments, and restore gravity to General Relativity. The code uses a multigrid relaxation scheme to solve the nonlinear field equation in the quasistatic approximation. We perform simulations of a flat selfaccelerating DGP model without cosmological constant. However, the type of nonlinear interactions of the branebending mode, which are the focus of this study, are generic to a wide class of braneworld cosmologies. The results of the DGP simulations are compared with standard gravity simulations assuming the same expansion history, and with DGP simulations using the linearized equation for the brane bending mode. This allows us to isolate the effects of the nonlinear selfcouplings of the field which are noticeable already on quasilinear scales. We present results on the matter power spectrum and the halo mass function, and discuss the behavior of the brane bending mode within cosmological structure formation. We find that, independently of CMB constraints, the selfaccelerating DGP model is strongly constrained by current weak lensing and cluster abundance measurements.
pacs:
95.30.Sf 95.36.+x 98.80.k 98.80.Jk 04.50.KdI Introduction
The observed acceleration of the universe Riess et al. (2007); Kowalski et al. (2008); Dunkley et al. (2008); Giannantonio et al. (2006); Pietrobon et al. (2006) poses a fascinating challenge to physicists and cosmologists: it points towards new physics at an unusually low energy scale (eV), or large length scale (Mpc). For this reason, no natural explanation appears to exist up to now. Many attempts have been made to go beyond the minimal explanation, a cosmological constant or vacuum energy. They can be broadly classified into two categories: acceleration is due to an additional, smooth stressenergy component with negative pressure (dark energy Frieman et al. (2008)); or it is caused by gravity itself, through modifications to General Relativity (GR) on cosmological scales. While there are strong constraints on deviations from GR in the Solar System Will (2006), gravity is remarkably weakly constrained on cosmological length scales. This provides an independent motivation for studying modified gravity models in the context of cosmology.
Smooth dark energy models have enough freedom to reproduce essentially any background expansion history of the universe. In order to distinguish modified gravity from the smooth dark energy scenario, it is thus necessary to study the growth of cosmological structure. A wealth of observables can be used for this purpose Zhang et al. (2007); Jain and Zhang (2007); Song and Koyama (2008); Knox et al. (2006); Schmidt (2008); Song (2005, 2006); Jain and Zhang (2007); Tsujikawa and Tatekawa (2008); Zhang (2006); Song et al. (2007); Schmidt et al. (2007), e.g. weak lensing, galaxyCMB crosscorrelation, galaxy cluster abundances, and many more. However, almost all of these measurements are affected by nonlinearities in the matter density field, or have most of their information on nonlinear scales. All viable modified gravity models include a nonlinear mechanism to restore General Relativity in highdensity environments, which is necessary in order to satisfy Solar System constraints. This mechanism has to be taken into account in order to make accurate predictions for observables on nonlinear scales in modified gravity. For one representative of modified gravity, Oyaizu et al. (2008); Schmidt et al. (2008) showed that the effects of this nonlinear chameleon mechanism are significant for models that satisfy Solar System constraints.
One of the most popular modified gravity scenarios is the DvaliGabadadzePorrati (DGP) model Dvali et al. (2000). In this model, matter and radiation are confined to a 4dimensional brane in a 5dimensional bulk spacetime. While gravity propagates in five dimensions on the largest scales, it becomes 4dimensional below a certain crossover scale . On scales smaller than , and when gravity is weak, DGP gravity can be described as an effective fourdimensional scalartensor theory Nicolis and Rattazzi (2004). The scalar field, referred to as brane bending mode, is massless but has nonlinear derivative interactions which suppress the field in highdensity environments, restoring GR locally and allowing DGP gravity to pass Solar System constraints.
Applied to homogeneous and isotropic cosmology, the DGP model allows for two solutions Deffayet (2001). In one branch of the theory, the selfaccelerating branch, the effective weakening of gravity on large scales leads to an accelerated latetime expansion of the brane, without any cosmological constant. The scalar branebending mode mediates a repulsive force, leading to a smaller effective gravitational constant on large scales. In the other solution, the normal branch, there is no accleration, and the branebending mode mediates an attractive force. In linear perturbation theory around its de Sitter limit, the brane bending mode in the selfaccelerating DGP model has been shown to have the wrong sign in the kinetic term (“ghost”) Luty et al. (2003); Nicolis and Rattazzi (2004), suggesting an instability, although the situation in the nonlinear case is not clear Dvali (2006); Koyama (2007). The normal branch is free of the ghost.
Linear perturbation theory around a cosmological background in DGP has been derived in, e.g. Koyama and Maartens (2006); Song et al. (2007); Cardoso et al. (2008). In this approach, it is possible to predict the expansion history, CMB anisotropies, and the ISW effect in DGP Fang et al. (2008a). Using this information, Fang et al. (2008b) have shown that the selfaccelerating DGP model without cosmological constant, which has only one free parameter, , is disfavored at the level by current CMB and Supernova data. This is due to the earlier onset of acceleration in DGP, and the additional suppression of growth through the branebending mode.
While the selfaccelerating branch of DGP is thus challenged by theory and observations, much work is being done on extending the DGP model to higher codimensions, in the context of degravitation Dvali et al. (2007); de Rham et al. (2008). These higherdimensional models are expected to bring the expansion history close to that of CDM, while exhibiting a similar effective scalartensor regime as in the original DGP model (see also Afshordi et al. (2008)). The form of the nonlinear interactions of the branebending mode are also generic to these generalized models (see also Nicolis et al. (2009)). Hence, it will be straightforward to perform simulations of these models once the expansion history and evolution of linear perturbations is worked out. In addition, the standard DGP model is able to satisfy the cosmological constraints if a cosmological constant (brane tension) is included Lombriser et al. (2009).
So far, studies of the DGP model in the context of cosmology have mostly dealt with the linearized theory, valid on large scales. However, as soon as perturbations in the matter density become of order unity, the nonlinear interactions of the brane bending mode become important (e.g., Koyama and Silva (2007)). Hence, in order to make predictions for observables in the nonlinear regime, the full nonlinear equations of motion must be solved in conjunction with the evolution of the matter perturbations. Recently, perturbative approaches have been developed to extend the predictions into the quasilinear regime Koyama et al. (2009). Also, Khoury and Wyman Khoury and Wyman (2009) have used a sphericallysymmetric approximation of the branebending mode interactions in their Nbody simulation of braneworld models. In this paper, we present the results of an Nbody simulation of the selfaccelerating DGP model, which selfconsistently solves the full nonlinear equation for the branebending mode and its effect on the motion of particles (we compare our results with the approximation used in Khoury and Wyman (2009) in Section VI.3). We solve the equation of motion in the quasistatic approximation, neglecting time derivatives with respect to spatial derivatives, as is usually done in Nbody simulations (Section IV.3 presents a consistency test of this approximation).
We present results on the matter power spectrum, the halo mass function (abundance of dark matter halos), and the behavior of the brane bending mode in the cosmological context. These results can be used to strengthen constraints on the selfaccelerating DGP model significantly, e.g., through observations of weak lensing and the abundance of galaxy clusters. More generally, the simulations presented here can serve as starting point for building a model of nonlinear structure formation in braneworld scenarios, and for benchmarking perturbative approaches Koyama et al. (2009).
In Section II, we describe the DGP model, and present the relevant equations and analytical test cases. Section III presents the code, which is benchmarked in Section IV. The cosmological simulations are described in Section V. Results and their impact on cosmological constraints on the DGP model are presented in Section VI and Section VII. We conclude in Section VIII.
Ii DGP model
ii.1 Background evolution
The DvaliGabadadzePorrati model Dvali et al. (2000) consists of a spatially threedimensional brane in a 4+1 dimensional Minkowski bulk. Matter and all interactions except gravity are confined to the brane. The gravitational action consists of the fivedimensional EinsteinHilbert action plus a term which leads to the 4D gravity limit on small scales:
(1) 
Here, , stand for the bulk coordinates and metric, while , are the induced coordinates and metric on the brane, and , denote the corresponding Ricci scalars. In the following, we will drop the notation where no confusion can arise. The boundary term is added to the action in order to ensure that variation with respect to leads to the correct fivedimensional Einstein equations (e.g., Gibbons and Hawking (1977)).
The two gravitational constants, or Planck masses , appearing in Eq. (1) can be related via a length scale, the crossover scale:
(2) 
On scales above , gravity becomes fivedimensional, with forces falling off as . Below , gravity is fourdimensional, but not Einsteinian gravity, a point to which we return below.
Since all matter is thought as confined to the brane, the fivedimensional metric has to obey nontrivial junction conditions over the brane Deffayet (2001). Assuming an empty Minkowski bulk, a spatially flat brane, and a homogeneous and isotropic matter distribution on the brane, the junction conditions lead to the following analogue of the Friedmann equation for the scale factor of the induced metric on the brane:
(3) 
The junction conditions leave two possible branches of the theory, determined by the sign on the left hand side of Eq. (3). In the following, we focus on the branch with the ’’ sign, which asymptotes to a latetime de SitterUniverse, , and is correspondingly called the accelerating branch. Then, in the matterdominated epoch, neglecting radiation, and again assuming no cosmological constant or curvature, Eq. (3) can be rewritten as:
(4)  
and is the average matter density today. This expansion history is clearly different from CDM, and corresponds to an effective dark energy with in the matterdominated era at high redshifts.
ii.2 Cosmological perturbations and brane bending mode
The propagation of light and particles on the DGP brane is completely determined by the perturbed 4D FriedmannRobertsonWalker metric:
(5) 
However, in order to determine the evolution of the metric potentials on the brane, it is necessary to solve the full 5D Einstein equations Koyama and Silva (2007); Song et al. (2007). An additional scalar degree of freedom associated with local displacements of the brane appears, the socalled branebending mode which couples to matter. In our convention, is dimensionless, instead of being scaled to the Planck mass .
In the decoupling limit of DGP Nicolis and Rattazzi (2004), when setting gravitational interactions to 0, the selfinteractions of the field remain constant. Hence, while perturbations of the metric higher than linear order can be neglected for cosmological studies, it is crucial to consider the selfinteractions in the branebending mode.
In the quasistatic regime, for scales , time derivatives can be neglected with respect to spatial derivatives, and the equation for the branebending mode reads (e.g., Koyama and Silva (2007)):
(6) 
where the derivatives are with respect to comoving coordinates , is the matter density perturbation, and the function (for the accelerating branch) is given by:
(7) 
Note that is always negative in the selfaccelerating branch. In Section IV.3, we show results of a consistency test of the quasistatic assumption used in Eq. (6).
The branebending mode influences the dynamics of particles through the dynamical potential , which, assuming the same boundary conditions for and , is given by:
(8) 
where the Newtonian potential satisfies the usual Poisson equation:
(9) 
In contrast, the propagation of photons, determined by the lensing potential , is not directly affected by .
ii.3 Weak brane regime
If the gradient of the field is small, i.e. for small gravitational accelerations, Eq. (6) can be linearized, yielding a standard Poisson equation:
(10) 
We will refer to a cosmology with this linearized equation as the linearized DGP model. In this regime, also called weakbrane phase, becomes proportional to the Newtonian potential , corresponding to a constant rescaling of the gravitational constant through Eq. (8):
(11) 
Substituting the linear solution into Eq. (6), we obtain a rough estimate for the overdensity at which the nonlinear interactions become important:
(12) 
For a selfaccelerating DGP model that leads to cosmic acceleration today, , , and the prefactor in Eq. (12) is of order unity today. Hence, the selfcoupling of the brane bending mode becomes important as soon as the matter density field becomes nonlinear, .
ii.4 Vainshtein effect
In general, Eq. (6) is difficult to solve in full generality due to the nonlinearity in the derivative terms. However, an instructive test case, a spherically symmetric matter distribution, can be solved analytically. In the spherically symmetric case, Eq. (6) becomes Koyama and Silva (2007):
(13)  
For simplicity, we assume a spherical mass of radius with uniform density, and set . Then, we can integrate Eq. (13) once and obtain the gravitational acceleration in DGP:
(14) 
where is the Newtonian acceleration of the spherical mass, and:
(15) 
Here, denotes a characteristic scale of the solution, the Vainshtein radius:
(16) 
and is the gravitational radius of the mass. For very large distances, , approaches the constant , which exactly matches the linear solution, Eq. (11). Note also that by substituting in Eq. (12), we see that the nonlinearity criterion is directly proportional to .
In the opposite limit, , becomes suppressed with respect to the linear solution, and approaches the small constant inside the mass (assuming ). This Vainshtein effect Vainshtein (1972); Deffayet et al. (2002) amounts to restoring GR in deep potential wells, while far away from the mass, gravity is in the scalartensor regime.
ii.5 Plane wave solution
Another exact solution to Eq. (6) exists: for a simple plane wave density perturbation, . In this case, the two nonlinear terms exactly cancel, and we are left with the linear solution, Eq. (10). In Appendix B this is used to test our code for perturbations of different wavelengths. Also, the plane wave and spherically symmetric solutions can be considered as two limiting cases for understanding the field behavior in the cosmological context. We will return to this point in Section VI.3.
Iii Code implementation
We use an extension of the code first described in Oyaizu (2008) for our simulations of DGP gravity. The code is a standard particlemesh Nbody code Hockney and Eastwood (1981); Klypin and Holtzman (1997) assuming collisionless dark matter only, with a fixed grid of size and periodic boundary conditions, and uses secondorder accurate leapfrog integration for the particle propagation. The major addition is a relaxation solver for the nonlinear equation of the brane bending mode , Eq. (6). The simulation proceeds by performing steps in the scale factor .
Except for the different expansion history, the equations for particle propagation are identical to those in standard CDM simulations, and are expressed in terms of comoving coordinates , :
(17)  
(18) 
Here, and are the position and momentum of particle in code units (see Appendix A), and is taken from Eq. (4). The main modification enters in the determination of . One time step of the Nbody simulation proceeds as follows:

The density field on a fixed grid is assigned using the particle positions at scale factor via the cloudincell method, i.e. each particle corresponds to a uniform density cube of side length and mass .

Using the same cloudincell interpolation as for particles, the gradient of is calculated at each particle position and is used to update the particle momenta from to .

Using the momenta at , the particle positions are updated from to , and the process starts from the beginning.
Thus, in each time step we have to solve for the branebending mode which contributes to the dynamical potential via Eq. (8). In other words, given the density field, we have to solve the nonlinear elliptical differential equation Eq. (6). We use a GaussSeidel relaxation scheme together with the Newton method. A crucial tool is the multigrid Brandt (1977, 1984); Briggs et al. (2000), a hierarchy of coarser grids which are essential to speeding up the convergence Oyaizu (2008): as the relaxation scheme, which operates locally, gets rid of smallscale errors very efficiently, the coarser grids quickly reduce the longwavelength error modes which are hard to eliminate on the fine grid alone. For details of the implementation, see Appendix A.
While the field is in general wellbehaved, the convergence properties become worse for strongly inhomogeneous density fields, because the nonlinearities in Eq. (6) are in the derivatives of the field. In order to reach the desired convergence in our cosmological simulations, it is necessary to smooth the density field entering the r.h.s. of Eq. (6) with a Gaussian kernel , with set to the size of a grid cell. This smoothes over the noise in the density field due to the discreteness of particles.
Iv Code tests
In this section, we present different tests our code was put through to benchmark its performance and limitations. We focus on tests applying to the modified gravity sector. For the results of standard Nbody code tests for this code, see Oyaizu (2008).
iv.1 Spherical mass test
In order to study how well the code reproduces the exact result for a spherical mass in DGP, we start out with a grid with grid cells and particles. We arbitrarily assume a box size of and set . All particles are moved into a spherical mass of radius and uniform density, which then corresponds to a mass of (assuming ). Given the density field assigned from the particle positions, we then use the relaxation solver to solve for , and measure the field values and radial acceleration throughout the box.
For the purposes of this test, we can vary the two parameters and in Eq. (6) independently. determines the strength of the additional force mediated by [e.g., Eq. (10)], which is attractive for positive and repulsive for . For given , controls the Vainshtein radius, i.e. the scale where the selfinteractions of the field become important.
Fig. 1 shows the field solution (left panel) and relative deviation of the acceleration from the Newtonian value (right panel) as a function of distance from the center, , for , grid cells, corresponding to , and . The agreement with the analytical solution (thick lines) is generally very good, except at large radii, where the periodic boundary conditions become important, and around , where artifacts of interpolating a sphere onto a cubic grid become visible. At distances of order the size of the grid cells, the force resolution becomes worse leading to a growing scatter in the measured acceleration. Note that the analytical solution necessarily assumes an isolated mass, since the superposition principle cannot be used due to the nonlinear field equation. For this reason, we added an arbitrary zeropoint to the field in order to match the analytical solution in the left panel of Fig. 1. Of course, such a zeropoint does not affect any observable quantity, such as the acceleration shown in the right panel.
We tested the field solution for several different values of , , , and . Given the caveats pointed out, we found very good convergence to the analytical solution. In particular, increasing reduces the interpolation artifacts at , while decreasing reduces the impact of the periodic boundary conditions. Solutions for different values of are shown in both panels of Fig. 1. While the solid lines show the exact solution to the full equation, Eq. (13), the dotted lines show the linearized result, Eq. (10). Clearly, for the parameters chosen, the effect of the nonlinearity is significant for of a few or less.
iv.2 Convergence and resolution tests
In most test cases we studied, the density field is sufficiently smooth that the relaxation converges to the desired low tolerance. However, for density fields with considerable smallscale inhomogeneities, such as the cosmological density field at late times, we find that this level is not achievable. This is because the nonlinearity of the equation is important as soon as the overdensity becomes of order unity, which in the Nbody simulation corresponds to particle per grid cell. Equivalently, the Vainshtein radius for a single particle in our Nbody simulation is of order the grid scale (see Fig. 2, bottom left). Thus, the field reacts much more strongly to noise due to the discreteness of particles than the Newtonian potential, leading to increased residual errors of the approximate solution.
If we write Eq. (6) in code units (see Appendix A) as , then for an approximate solution , is the dimensionless residual, where is the size of a grid cell. In the following, we use the RMS of the dimensionless residual, as a benchmark for the convergence of the field solution. We performed tests with sine wave density fields of various wavelengths in order to determine what residual is acceptable in our simulations (see Appendix B). In case of a pure sine wave, the nonlinearity in Eq. (6) vanishes, and the exact nonlinear solution is identical to the linearized solution. We found that residuals of are safe, as for residuals at that level, the errors in the solution are negligible compared to the unavoidable truncation errors from taking numerical derivatives on the grid.
In order to reduce the noise in the density field in the cosmological simulations and improve the solution to the acceptable level of , we increase the number of particles from to , and smooth the density field entering the r.h.s. of Eq. (6) as described in Section III. With these steps, the solution converges with a residual for all box sizes (left panel of Fig. 2). Note that the dimensionless RMS residuals never exceed a fraction of of the RMS of the field solution.
Since the force resolution in our simulations is in any case limited to scales above , a smoothing on the grid scale is not expected to degrade the resolution of the simulations significantly. We checked this by running linearized DGP simulations given by Eq. (10) (Section II.3) using the same smoothing of the r.h.s. of Eq. (6), and comparing the resulting power spectrum to that of linearized DGP simulations without smoothing, for the same initial conditions. The result is shown in the right panel of Fig. 2. As expected, the effect of smoothing increases towards larger . The smoothing effect on the matter power spectrum is positive, since the smoothing removes power in the field on the smallest scales, and mediates a repulsive force. This leads to an enhancement in the matter power spectrum.
The smoothing effect remains below 4% for , which we adopt as the maximum wave number considered for each box, where is the Nyquist frequency of the grid. From our studies with different smoothing radii , we found that the smoothing effects on full DGP simulations show a similar dependence as those of the linearized DGP simulations, but are smaller by a factor of . This is understandable since the field is suppressed in dense regions in the full simulations, reducing the effect of the smoothing. Taking into account this factor, we correct the power spectra measured in the full DGP simulations for the smoothing effects using the curves shown in Fig. 2 (right panel). Note that in any case these effects are at the level of few percent or less.
In order to assess the effect of the finite grid resolution, we also performed simulations with and , and . Fig. 2 (right panel) also shows the deviation of the power spectrum in these lowresolution simulations from the simulation with the same initial conditions. The vertical lines indicate for each case, i.e. the maximum we consider for each simulation box. Below , the deviations are less than 10% for either or . Note also that the resolution effects are independent of the type of simulation (DGP, linearized DGP, or GR). Hence, they cancel when measuring the deviation of between different simulation types, and we expect this deviation to be accurate to the percent level for the range in we consider.
iv.3 Quasistatic approximation
In solving for the brane bending mode , we have assumed that all timederivative terms in the full equation of motion of (in terms of physical coordinates):
(19) 
can be neglected with respect to the spatial derivatives in Eq. (6). Up to second order, there are three time derivative terms which have been neglected:
(20)  
(21)  
(22) 
and appear as in the equations of motion. While we cannot rigorously prove that the quasistatic approximation holds, we can perform a consistency test by measuring the terms Eqs. (20)–(22) in the simulations, and checking whether they are indeed small. The time derivatives of are calculated in each grid cell using the two neighboring time steps. We then take the RMS of each term over the grid, , and compare with the RMS of the spatial Laplacian, . Fig. 3 presents and relative to the spatial derivatives for our largest and smallest box sizes, and shows that they are well under control. The relative magnitude of all time derivatives are of order a few or less. This is in keeping with the expectation that , and in our simulations. We conclude that the quasistatic approach is a selfconsistent approximation.
V Cosmological simulations
0.258  
0.138  
66.0  
2.37  
0.0888  
0.0954  
0.998  
()  
0.6566 
[]  
# of  QCDM  6  6  6  6 
boxes  Lin. DGP  6  6  6  6 
Full DGP  6  6  6  6  
[]  0.50  0.79  1.57  3.14  
[]  0.78  0.50  0.25  0.13  
[]  219  57.3  7.17  0.90 
We performed a suite of cosmological simulations for three different types of gravity: unmodified GR, corresponding to a smooth dark energy model with the same expansion history as DGP, referred to as QCDM; linearized DGP gravity, using the linearized equation of motion [Eq. (10)], corresponding to a timedependent rescaling of Newton’s constant; and full DGP, solving for the full nonlinear solution [Eq. (6)]. Comparing linearized DGP with full DGP allows us to study the Vainshtein effect in a cosmological setting.
The cosmological parameters are those of the bestfit flat selfaccelerating DGP model to WMAP 5yr data Fang et al. (2008c) and are summarized in Tab. 2. We generated the initial conditions at () from a modified transfer function output of CAMB Lewis et al. (2000) for a flat CDM model with . The CDM transfer function was corrected for small earlytime modified gravity effects in the DGP model using the PPF approach, as detailed in Appendix C.
The simulations were run with grid cells on a side, and particles, i.e. one particle per grid cell, to reduce the shot noise in the density field (see Section IV.2). We performed six simulations each for four different box sizes, from up to (Tab. 2). On an 8core machine, the QCDM and linearized DGP runs require 10h of computing time, while the full DGP simulations require about 300h.
Fig. 4 (left panel) shows the combined power spectrum from all box sizes for the QCDM simulations, including bootstrap errors. Also shown is the nonlinear power spectrum for QCDM, calculated from the linear power spectrum using the halofit procedure Smith et al. (2003). The power spectrum measured in each box is used up to (see Tab. 2), and different boxes are combined weighting by volume. The lower panel of Fig. 4 (left) shows the power spectrum relative to the halofit prediction, measured separately for each box size. The power spectra measured in different boxes clearly agree within the errors, and the deviations from the halofit prediction are within the accuracy (%) expected from this fitting procedure (especially given the significant departures from CDM of the simulated expansion history).
Fig. 4 also shows the nonlinear power spectrum for a flat CDM cosmology fixed to the same , , and primordial normalization, corresponding to a linear normalization of today. Note that due to the different expansion history and the earlier onset of acceleration in DGP, the linear growth factor is suppressed by % in QCDM relative to CDM. The repulsive branebending mode in DGP suppresses growth further, leading to a suppression in the linear regime of % relative to CDM. In the following, we will always compare the DGP results with QCDM, so that the expansion history effects are taken out, and all deviations are strictly due to modifications of gravity.
We have also measured the mass function of dark matter halos in our simulations. The spherical overdensity halofinding and combination of different simulation boxes was done as described in Schmidt et al. (2008). As our grid resolution is the same as in the simulations analyzed in Schmidt et al. (2008), while we have a factor of 8 more particles, we conservatively increase the minimum required number of particles from 800 to 6400. The corresponding mass thresholds are given in Tab. 2. We define the halo mass as the mass enclosed within a sphere of radius , so that the average density within the sphere is 200 times the average matter density in the universe, .
Fig. 4 (right panel) shows the mass function, , measured in QCDM simulations, in comparison with the prediction calculated from the linear QCDM power spectrum using the fitting formula from Tinker et al.Tinker et al. (2008). The lower panel shows the relative deviations from the prediction, separately for each simulation box. Within the mass range accessible with our simulations, , the agreement with the prediction and among different boxes is good. We also show the predicted mass function for a CDM model with the same initial power spectrum in Fig. 4. Apparently, the number of halos above is significantly reduced due to the suppressed growth in QCDM.
Vi Results
In order to get a visual impression of some of the physics in DGP Nbody simulations, we show slices through one simulation box at () in Fig. 5. The slices are 64 cells thick, and for each pixel we take the maximum absolute values of each quantity over the thickness of the slice, for better visibility. The density field (top left; in logarithmic scale) is difficult to distinguish visually from the QCDM result for the same run, as the effects on the power spectrum are at the % level (Section VI.1). The top right panel in Fig. 5 shows the dynamical potential , exhibiting the potential wells of massive collapsed structures. The branebending mode is shown in the lower left panel. On linear scales, is proportional to the potential (Section II.3), but evidently does not follow the potential within deeper potential wells, making it appear smoother. To make this more clear, we show the difference of the solution from the linear value [Eq. (10)], , in the lower right panel of Fig. 5: is suppressed () in overdense regions, showing the Vainshtein effect at work in a cosmological setting. Note that quite lowmass structures which are not conspicuous in the potential already lead to a suppression of . See Section VI.3 for a discussion of the Vainshtein effect in dark matter halos.
vi.1 Power spectrum
Fig. 6 (left panel) shows the relative deviation of the linearized and full DGP power spectra from the QCDM result at redshift 0. The error bars are obtained from the 6 separate runs using a bootstrap procedure. By comparing simulation runs with the same initial conditions, and then averaging the deviations, most of the cosmic variance cancels out and we are able to obtain significantly less scatter. Also shown is the predicted deviation in the linear power spectrum. Note first that in the selfaccelerated DGP branch, the scalar field mediates a repulsive force, leading to a suppression of the growth of structure. Due to the scaleinvariant modification of gravity in linearized quasistatic DGP [Eq. (10)], the predicted linear deviation is also scaleindependent. On linear scales, , both linearized and full DGP simulations agree well with the linear prediction.
Apparently, the full DGP result departs significantly from linearized DGP on quasilinear to nonlinear scales. The Vainshtein effect begins to operate wherever overdensities become of order unity, and suppresses the deviation from GR. Note that small effects of the nonlinear equation can already be seen on quite large scales, corresponding to , not far from the acoustic features in the matter power spectrum. While we do not expect dramatic effects on cosmological parameter constraints from BAO, a modeling of these nonlinear effects will be necessary for precision constraints in the context of DGP and similar braneworld models.
We also show the deviation of from , where is calculated from the corresponding linear power spectrum using the standard halofit prescription. halofit describes the linearized DGP power spectrum reasonably well up to , in agreement with the findings of Laszlo and Bean (2008), while it follows neither the linearized nor full DGP at higher . We have not explored whether phenomenological modifications of halofit allow for an improvement of the fit.
The dashed points at high extend the range in frequency up to of our smallest box, . While resolution effects are still expected to cancel out to first order in this representation, these points only serve to indicate that the trends seen for full and linear DGP at lower continue towards smaller scales.
The right panel of Fig. 6 shows the evolution of the modified gravity effects on the power spectrum as function of redshift. On large scales , the deviation is almost scalefree and evolves as predicted by linear theory. At earlier times, the density field is nonlinear only on smaller scales. Hence, the Vainshtein effect becomes visible in the power spectrum only at higher . However, it does affect the power spectrum deviation significantly already at .
vi.2 Halo mass function
The abundance of massive dark matter halos is a sensitive probe of the growth of structure in the Universe White et al. (1993); Eke et al. (1998); Borgani et al. (2001); Vikhlinin et al. (2009). In particular, the number of the most massive halos which host galaxy clusters depends exponentially on the amplitude of the matter power spectrum. Fig. 7 shows the effect of the modification of gravity in DGP on the halo mass function, relative to the QCDM effective dark energy cosmology. We investigated the effect of the smoothing used in the full DGP simulations by comparing linearized DGP simulations with and without smoothing, as done for the power spectrum (Section IV.2). Above our rather conservative mass threshold for each box, we only found a small effect of the smoothing on halo masses in the linearized DGP simulations. Since the field is in any case suppressed within halos due to the Vainshtein mechanism, we expect the smoothing effects to be even smaller for the full DGP simulations.
As expected, the full DGP simulations are somewhat closer to QCDM than the linearized DGP case in Fig. 7. The suppression in the mass function is significant for , reaching more than 30% for massive clustersize halos. The suppression is smaller for lowermass halos. This is presumably because these halos formed earlier in cosmic history, when the modified gravity effects of DGP were significantly smaller.
Note that when compared to CDM, this suppression comes in addition to the larger suppression of the mass function due to the expansion history of DGP (see dashed line in Fig. 7). However, we expect the effect of the DGP modification of gravity to be fairly independent of the expansion history. In particular, we expect a similar effect on the halo mass function in generalized braneworld models which exhibit an expansion history close to CDM de Rham et al. (2008); Khoury and Wyman (2009). Hence, the mass function is expected to be a sensitive probe for braneworld modified gravity models, as is the case for gravity Schmidt et al. (2008). Note that for the normal branch of braneworld gravity, the sign of the effect in Fig. 7 will be reversed, leading to an enhancement of the number of massive halos.
vi.3 Brane bending mode and Vainshtein effect
Fig. 8 (left panel) shows the average profiles of the branebending mode , the dynamical potential , and the lensing potential of dark matter halos in the full DGP simulations. The branebending mode has been scaled by , so that it agrees with for the linearized solution. We only use our highest resolution simulation box for this measurement (). The inner radius of the profiles is given by one grid cell. The profiles were measured by spherically averaging many lines of sight for each halo with , and then stacked, scaling each profile to the radius of the respective halo. This was done in order to reduce the significant scatter in the profiles.
While the lensing potential obeys the standard Poisson equation in DGP, the dynamical potential receives a contribution from the brane bending mode [Eq. (8)]. Apparently, the field flattens out towards the halo center, which is in qualitative agreement with the nonlinear field solution for a spherical mass (see, e.g. Fig. 1). Note that turns away from the linear solution when the overdensity is of order a few, in agreement with the estimate in Section II.3.
The quantity which is actually observable however is the gradient of the field, shown in Fig. 8 (right panel), which gives essentially the deviation of the acceleration from the Newtonian acceleration. For each halo, we measure the radial gradient relative to the halo center. We again scaled by to match the linearized solution to . From the halo exterior towards the interior, the gradient of grows initially, turns around at , and shrinks for even smaller radii. As is clearly shown by comparing the dynamical potential with its Newtonian value, , any observable deviations from GR are indeed suppressed within massive halos. We find that for halos between and , the deviation in the acceleration is around in the inner parts, while it is suppressed down to for halos around . This trend with mass is expected, given that for the spherically symmetric solution.
In order to model the behavior of the brane bending mode in the cosmological context, we can make use of the two limiting cases presented in Section II. One is the spherically symmetric solution (Section II.4), which can be used in modeling a spherical collapse of dark matter halos in DGP Lue et al. (2004). Moreover, Khoury and Wyman Khoury and Wyman (2009) have based the approximate field solution in their Nbody simulation on this limiting case. More realistically however, cosmic structure forms by collapsing into nonradially symmetric structures such as filaments. The plane wave density perturbation can serve as another, albeit rather extreme, test case. For such a perturbation, the nonlinear selfcoupling of vanishes (Section II.5). Hence, we expect the Vainshtein effect to be typically weakened in a nonradially symmetric setting, compared to the spherically symmetric case.
In their ansatz, Khoury and Wyman assumed that the spherically symmetric solution Eqs. (14)–(15) holds wherever the field becomes nonlinear. Since in this solution, the gravitational constant is rescaled locally by , the Poisson equation for can be written in this ansatz as:
(23) 
where is now a function of the local overdensity. Substituting in Eq. (16) [or Eq. (12)] we see again that , with an order unity coefficient. Khoury and Wyman Khoury and Wyman (2009) chose:
(24) 
Hence, in the DGP model we simulated, . Given the density field in our DGP simulations, we can solve Eq. (23) by Fourier transform as done in Khoury and Wyman (2009), subtract the Newtonian potential, and compare the brane bending mode from this ansatz with our numerical solution of the full equation.
Fig. 8 (right panel) shows the result for stacked halo profiles of the radial gradient of . While shows a qualitatively similar behavior to the full solution, the suppression due to the Vainshtein effect appears to set in at significantly larger radii in than in the full solution. This is qualitatively in agreement with our expectation that nonradially symmetric configurations experience a weaker nonlinear suppression. The fact that does not approach the linear solution at large radii is a result of the approximation Eq. (23), as was derived in the appendix of Khoury and Wyman (2009). The precise tensorial structure of Eq. (6) restores to the linear solution at large distances . It would be worth investigating whether this simplified (and computationally much less expensive) approach can be extended to recover the linear solution at large scales, an essential feature of the full branebending mode interactions.
Although we did not compare our results with a full simulation based on the spherically symmetric approach of Khoury and Wyman (2009), these results seem to indicate that this approximation overestimates the nonlinear suppression of the branebending mode, which might affect observables such as the power spectrum or halo mass function. We point out however that the crossover scale in the braneworldinspired model considered in Khoury and Wyman (2009) is an order of magnitude smaller than our , and the nonlinearity in Eqs. (23)–(24) only becomes effective at much higher overdensities , so that for the models studied in Khoury and Wyman (2009), the solution is possibly less sensitive to this approximation.
Vii Constraints on the selfaccelerating DGP model
We now briefly discuss the impact of our simulation results on cosmological constraints on the selfaccelerating DGP model without . First, regarding, baryon acoustic oscillations, our results for the DGP power spectrum on quasilinear scales show that nonlinear effects are not significantly enhanced at the BAO scale. We estimate that any modified gravity corrections to the BAO scale are below the percent level, and hence well within the uncertainties of current BAO measurements Eisenstein et al. (2005); Percival et al. (2007). Including baryon acoustic oscillations will increase the power of the combined CMB and Supernova constraints on the selfaccelerating DGP model Fang et al. (2008b) with nonzero curvature to Hu ().
Second, weak lensing measurements constrain the amplitude of the matter power spectrum today, which in this model is significantly smaller due to the suppression of growth by the earlier onset of acceleration and the repulsive force mediated by the branebending mode. The linear power spectrum normalization at in the models we simulated is and , while for a CDM model with the same primordial normalization, it is . We found that the nonlinear matter power spectrum in the full DGP simulations is always below that of QCDM (up to , Fig. 6). Hence, the suppression of the nonlinear power spectrum in the selfaccelerating DGP model corresponds to a reduction in the inferred linear normalization today of at least 0.1, with respect to a CDM model with the same initial conditions. Such a deviation is disfavored by about 1.5 standard deviations by current weak lensing measurements Fu et al. (2008).
The abundance of massive dark matter halos is probed by cluster surveys. We showed in Section VI.2 that the abundance of halos above is suppressed by around 20% relative to QCDM, which itself only predicts half as many halos in that mass range as a CDM model with the same primordial power. The latter again corresponds to a shift in by 0.1, which is constrained to more than by Xray cluster measurements Vikhlinin et al. (2009), when taking into account the systematic errors.
Viii Conclusions
The Nbody simulations of DGP gravity we presented here show how the selfinteractions of the brane bending mode , which are responsible for restoring General Relativity in dense environments, influence the formation of structure in the universe. In the selfaccelerating DGP model we simulated, this scalar field mediates a repulsive force, effectively weakening gravity. We indeed find that the field solution turns away from the linearized solution whenever the matter overdensity becomes of order a few, and that the repulsive force is suppressed by more than an order of magnitude compared to its linearized value in the center of massive halos. We also compared our full solution to that obtained in the approximate ansatz of Khoury and Wyman (2009). While their ansatz agrees well with our results in the densest regions, it seems to overestimate the nonlinear suppression in less dense regions, such as the outer regions and environments of halos. This is in line with our finding that the nonlinear suppression of is weaker for nonspherically symmetric configurations.
The nonlinear matter power spectrum in the DGP simulations shows that the suppression due to the repulsive field is amplified on nonlinear scales for the linearized DGP simulations. In the full DGP simulations, this suppression is smaller, and eventually turns around on Mpc scales. The deviations between linerized and full DGP power spectra are noticeable already on quasilinear scales, . At the BAO scale, modified gravity effects on the power spectrum are at the percentlevel.
We found that the abundance of massive dark matter halos is significantly suppressed in the DGP simulations, compared to a standard gravity simulation with the same expansion history. Again, the nonlinear suppression of the field alleviates this suppression, and hence has to be taken into account when using observations to constrain this type of modified gravity. The effect on the halo mass function is in fact large enough to make this an interesting observational probe of more general braneworld scenarios.
Independently of the CMB constraints Fang et al. (2008b), our results on the nonlinear structure formation strongly constrain the selfaccelerating DGP model (without ). In the future, we plan to extend our simulations to the normal branch of DGP, and more general braneworld(inspired) models de Rham et al. (2008); Khoury and Wyman (2009); Nicolis et al. (2009). The key part of the simulations, solving for the nonlinear interactions of the branebending mode, is generic to all of these models, and hence the code will be readily generalized to these cases. We also plan to attempt a modeling of the modified gravity effects, in the context of the halo model for example. This model can then serve as a framework for cosmological constraints on braneworld gravity which take into account the nonlinear mechanisms inherent to this model consistently. Simulations of both gravity and DGP have shown that most interesting phenomena appear on scales of 1 to tens of Mpc, and in the abundance and environments of clusters. Fortunately, these scales are well accessible to observations, e.g., through weak lensing, the Lyman forest, and cluster surveys, which can then be used as precision probes of gravity on cosmological scales.
Acknowledgements.
I am indebted to Scott Dodelson, Wayne Hu, and Andrey Kravtsov for invaluable input and guidance. I would like to thank Hiro Oyaizu and Marcos Lima for our previous collaborative work on the code and analysis, and Angela Olinto, Mike Gladders, Cora Dvorkin, Sam Leitner, Michael Mortonson, and Amol Upadhye for discussions. The support of the Fermilab computing staff is gratefully acknowledged. The simulations used in this work have been performed on the Joint Fermilab  KICP Supercomputing Cluster, supported by grants from Fermilab, Kavli Institute for Cosmological Physics, and the University of Chicago. This work was supported by the Kavli Institute for Cosmological Physics at the University of Chicago through grants NSF PHY0114422 and NSF PHY0551142.Appendix A Code implementation and discretization
This appendix describes details of the Nbody code implementation, focusing on the relaxation solver for Eq. (6) as the main nonstandard part of the code. The code is written in C++, uses OpenMP for parallelization of most timecritical operations, and employs FFTW fft () for the Fast Fourier Transforms.
a.1 Code units
The comoving code units follow the convention used in Shandarin (1980); Kravtsov et al. (1997), where
(25) 
and
(26) 
In the above definitions, is the comoving simulation box size in , is the number of grid cells in each direction, is the Hubble parameter today, is the critical density today, and is the fraction of nonrelativistic matter today relative to the critical density.
Transformed to code units, the equation for becomes:
(27) 
Here, acts with respect to code coordinates, , , and is the matter overdensity on the grid, determined from the particle positions. Note that we have not yet scaled to . After solving Eq. (27), is added to the standard Newtonian potential, which is solved for using a Fast Fourier Transform, to obtain the dynamical potential in DGP:
(28) 
a.2 Discretization of equation
We employ standard secondorder symmetric differences for the discretization of the second derivatives in Eq. (27):
(29)  
(30) 
and correspondingly for derivatives with respect to , . Here, stand for the grid indices, and the step size for the base grid and for the grid of refinement level . The finite differences are evaluated before squaring in calculating the nonlinear terms in Eq. (27).
We have tried different discretizations, such as going to higher order in the finite differences, and solving for the deviation of from the solution of the linearized equation [Eq. (10)], instead of solving for itself. The performance of those discretizations is comparable to or worse than the simple discretization Eqs. (29)–(30). This is understandable, since going to higher order amounts to making the derivative operations less local in order to be sensitive to error modes of longer wavelength. However, in our multigrid relaxation scheme, this is already done efficiently by the coarser grids, so going to higher order in a single relaxation does not improve performance. Hence, we stay with the straightforward discretization, Eqs. (29)–(30).
a.3 Relaxation algorithm
In general, a relaxation scheme operates by iteratively obtaining a better approximation to the solution given a previous guess, . On the grid, the solution is updated successively for each cell :
(31) 
where is the discretized field equation solved for . In case of a linear field equation, can be solved for straightforwardly. However, in our case the field equation is nonlinear in , and we instead solve for iteratively via Newton’s method. Writing the field equation Eq. (27) as , we need to solve:
(32) 
for , where we have suppressed the dependences of on neighboring grid cells. By expanding in a Taylor series around the current approximation, one step of Newton’s method works by updating with:
(33) 
Our relaxation step is thus given by Eq. (33), where is determined from Eq. (27) and the discretization Eqs. (29)–(30). In practice, the order in which we loop over grid cells in performing the relaxation Eq. (33) is important, since each relaxation step depends on the neighboring grid cells. In particular, we need to take care when parallelizing the relaxation. We use a generalized redblack scheme (see, e.g., Briggs et al. (2000); Oyaizu (2008)), by successively running over cells with modulo 4 = , . This was done in order to break dependences due to neighboring cells within each set, in particular due to the mixed derivatives Eq. (30), so that each set can be efficiently parallelized. We experimented with different ordering schemes and found that this choice performed best.
Now, assume we have an approximate solution to the field equation Eq. (27) on the fine grid , which differs from the true solution by the error , . Further, the residual is defined as . Since by assumption, we obtain the following residual equation:
(34) 
We expect that the relaxation on the fine grid has removed most smallscale error modes, so mainly consists of longer wavelength modes. Thus, we will solve for on the coarser grid , and then correct the approximate solution for this error. On the coarse grid , Eq. (34) reads:
(35) 
Here, the superscripts denote which grid a given quantity is defined on, and is the restriction operator mapping fields from the fine grid to the coarse grid. Note that Eq. (35) is of the same form as the equation on the coarse grid, , except with a different righthand side. Hence, Eq. (35) is solved for using the same algorithm as the original equation, but at one grid level higher. Once Eq. (35) is solved, we can correct for the error:
(36) 
where is the interpolation operator mapping the error from the coarse grid to the fine grid.
In summary, the multigrid relaxation proceeds as follows:
Thus, the relaxation proceeds as a nested loop going down to the coarsest grid, and then correcting the solutions on the successively finer grids. One such iteration through all grid levels is called a Vcycle. The coarsest grid we use has 4 cells on a side, which corresponds to a refinement of for a base grid.
For the interpolation operator , we use standard bilinear interpolation (consistent with the cloudincell scheme used for assigning densities and measuring accelerations on the grid). For the restriction , we use fullweighting, the transpose of the bilinear interpolation operator. See the appendix in Oyaizu (2008) for an explicit definition.
In order to ensure convergence at all times, we adopt a large value of , at the expense of computing time. Usually, one Vcycle reduces the residual by 1–2 orders of magnitude. We stop the relaxation when either is reached, or the convergence stalls at a certain level of . Either of this usually happens within 4 Vcycles. The value of for each time step is logged, allowing for a monitoring of the convergence status in the simulations.
Appendix B Sine wave test
In this section, we use the fact that the full solution is identical to the linearized solution for a plane wave perturbation (Section II.5) to test our code for perturbations on different length scales. We consider a simple sine wave density perturbation, with vector chosen to lie along the axis:
(37) 
In this case, the solution to the nonlinear equation is identical to the linear solution:
(38) 
Hence, the exact acceleration is given by:
(39) 
For a given level of residuals , we now demand that the errors due to the approximate solution of the equation of motion are small compared to the unavoidable truncation errors which are encurred by taking the gradient of the potential on the grid.
In the following, we set , , and the number of grid cells is set to . The truncation errors are measured by assigning the analytical solution for the potential to each grid point, and comparing the exact acceleration to the one obtained by derivation and interpolation from the grid. The truncation errors for two sine waves of different wavelength are shown as deviations in the acceleration in Fig. 9 (red points). We scale the deviations to the RMS acceleration of the sine wave, , in order to avoid far outliers obtained near the zeros of the sine wave when scaling to the exact solution at each point. For shorter wavelength modes, the truncation error increases and reaches order unity for waves on the Nyquist scale, as expected.
In order to estimate the additional error in the acceleration made due to insufficient convergence of the solution, we compare the acceleration measured from the nonlinear field solution coming out of the relaxation solver, to the corresponding solution of the linearized equation Eq. (10). We vary the parameters and which control the strength of the nonlinearity in Eq. (6). In Fig. 9, we show the measured deviations in the acceleration from the linear solution for different values of and , and in particular for the largest values that reached convergence. The final residuals of the solution in each case are also shown.
Apparently, the errors in the acceleration are completely negligible for residuals , consistent with results we found for other test cases such as the spherical mass. For sufficiently strong nonlinearity and short wavelengths, residuals of order are reached which lead to measurable deviations in the acceleration. However, for the residual values up to that we explored, the error due to the incomplete convergence is very small compared to the truncation error in all cases. Hence, we conclude that for residuals smaller than the conservative bound of , the solution we obtain is sufficiently accurate in order not to bias the particle dynamics.
Appendix C Correcting Initial Conditions for DGP
Due to the additional term in the Friedmann equation, the selfaccelerating DGP model starts to deviate from a CDM expansion history at relatively high redshifts. In addition, there are modifications to the growth of horizonscale modes. For this reason, even at the initial redshift of our simulations, , there are small departures in the matter transfer function in DGP as compared to a CDM model with the same earlyUniverse parameters .
In order to take these differences into account, we calculate the dark matter transfer function at for DGP and QCDM, using the PPF approach described in Hu and Sawicki (2007); Fang et al. (2008b), and for the corresponding flat CDM model. We do not include the effects of radiation on the transfer function, which is not necessary since we are only interested in the relative deviation of the DGP from CDM. Fig. 10 shows the relative deviation in from CDM for DGP and QCDM. Both QCDM and DGP transfer functions are slightly suppressed on all scales due to the effects of the earlier onset of acceleration. On superhorizon scales, the transfer function is further suppressed in DGP caused by the transition to 5D gravity on very large scales. In contrast, is less suppressed in QCDM because of the fluctuations in the effective dark energy.
We correct the transfer function obtained from CAMB for these small differences at by multiplying with , where is a simple function fit to the DGP curve shown in Fig. 10. Hence, these initial conditions are not quite correct for the QCDM simulations. However, the differences are very small on the scales probed by our simulations ().
Footnotes
 Linear power spectrum normalization today of a CDM model with the same primordial normalization.
References
 A. G. Riess, L.G. Strolger, S. Casertano, H. C. Ferguson, B. Mobasher, B. Gold, P. J. Challis, A. V. Filippenko, S. Jha, W. Li, et al., Astrophys. J. 659, 98 (2007), eprint arXiv:astroph/0611572.
 M. Kowalski, D. Rubin, G. Aldering, R. J. Agostinho, A. Amadon, R. Amanullah, C. Balland, K. Barbary, G. Blanc, P. J. Challis, et al., ArXiv eprints 804 (2008), eprint 0804.4142.
 J. Dunkley, E. Komatsu, M. R. Nolta, D. N. Spergel, D. Larson, G. Hinshaw, L. Page, C. L. Bennett, B. Gold, N. Jarosik, et al., ArXiv eprints 803 (2008), eprint 0803.0586.
 T. Giannantonio, R. G. Crittenden, R. C. Nichol, R. Scranton, G. T. Richards, A. D. Myers, R. J. Brunner, A. G. Gray, A. J. Connolly, and D. P. Schneider, Phys. Rev. D 74, 063520 (2006), eprint arXiv:astroph/0607572.
 D. Pietrobon, A. Balbi, and D. Marinucci, Phys. Rev. D 74, 043524 (2006), eprint arXiv:astroph/0606475.
 J. A. Frieman, M. S. Turner, and D. Huterer, Ann. Rev. Astron. Astroph. 46, 385 (2008), eprint 0803.0982.
 C. M. Will, Living Reviews in Relativity 9, 3 (2006), eprint arXiv:grqc/0510072.
 P. Zhang, M. Liguori, R. Bean, and S. Dodelson (2007), eprint arXiv:0704.1932 [astroph].
 B. Jain and P. Zhang, ArXiv eprints 709 (2007), eprint 0709.2375.
 Y.S. Song and K. Koyama, ArXiv eprints 802 (2008), eprint 0802.3897.
 L. Knox, Y.S. Song, and J. A. Tyson, Phys. Rev. D 74, 023512 (2006).
 F. Schmidt, Phys. Rev. D 78, 043002 (2008), eprint 0805.4812.
 Y.S. Song, Phys. Rev. D 71, 024026 (2005), eprint arXiv:astroph/0407489.
 Y.S. Song, ArXiv Astrophysics eprints (2006), eprint astroph/0602598.
 S. Tsujikawa and T. Tatekawa, ArXiv eprints 804 (2008), eprint 0804.4343.
 P. Zhang, Phys. Rev. D 73, 123504 (2006), eprint arXiv:astroph/0511218.
 Y.S. Song, I. Sawicki, and W. Hu, Phys. Rev. D 75, 064003 (2007).
 F. Schmidt, M. Liguori, and S. Dodelson, Phys. Rev. D 76, 083518 (2007), eprint arXiv:0706.1775.
 H. Oyaizu, M. Lima, and W. Hu, Phys. Rev. D 78, 123524 (2008), eprint 0807.2462.
 F. Schmidt, M. Lima, H. Oyaizu, and W. Hu, ArXiv eprints (2008), eprint 0812.0545.
 G. Dvali, G. Gabadadze, and M. Porrati, Physics Letters B 485, 208 (2000), eprint arXiv:hepth/0005016.
 A. Nicolis and R. Rattazzi, Journal of High Energy Physics 6, 59 (2004), eprint arXiv:hepth/0404159.
 C. Deffayet, Physics Letters B 502, 199 (2001), eprint arXiv:hepth/0010186.
 M. A. Luty, M. Porrati, and R. Rattazzi, Journal of High Energy Physics 2003, 029 (2003), URL http://stacks.iop.org/11266708/2003/i=09/a=029.
 G. Dvali, New Journal of Physics 8, 326 (2006), eprint arXiv:hepth/0610013.
 K. Koyama, ArXiv eprints (2007), eprint 0709.2399.
 K. Koyama and R. Maartens, Journal of Cosmology and AstroParticle Physics 1, 16 (2006), eprint arXiv:astroph/0511634.
 A. Cardoso, K. Koyama, S. S. Seahra, and F. P. Silva, Physical Review D (Particles, Fields, Gravitation, and Cosmology) 77, 083512 (pages 16) (2008), URL http://link.aps.org/abstract/PRD/v77/e083512.
 W. Fang, W. Hu, and A. Lewis, Phys. Rev. D 78, 087303 (2008a), eprint 0808.3125.
 W. Fang, S. Wang, W. Hu, Z. Haiman, L. Hui, and M. May, Phys. Rev. D 78, 103509 (2008b), eprint 0808.2208.
 G. Dvali, S. Hofmann, and J. Khoury, Phys. Rev. D 76, 084006 (2007), eprint arXiv:hepth/0703027.
 C. de Rham, S. Hofmann, J. Khoury, and A. J. Tolley, Journal of Cosmology and AstroParticle Physics 2, 11 (2008), eprint 0712.2821.
 N. Afshordi, G. Geshnizjani, and J. Khoury, ArXiv eprints (2008), eprint 0812.2244.
 A. Nicolis, R. Rattazzi, and E. Trincherini, Phys. Rev. D 79, 064036 (2009), eprint 0811.2197.
 L. Lombriser, W. Hu, W. Fang, and U. Seljak, ArXiv eprints (2009), eprint 0905.1112.
 K. Koyama and F. P. Silva, Phys. Rev. D 75, 084040 (2007), eprint arXiv:hepth/0702169.
 K. Koyama, A. Taruya, and T. Hiramatsu, ArXiv eprints (2009), eprint 0902.0618.
 J. Khoury and M. Wyman, ArXiv eprints (2009), eprint 0903.1292.
 G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2752 (1977).
 A. I. Vainshtein, Physics Letters B 39, 393 (1972), ISSN 03702693.
 C. Deffayet, G. Dvali, G. Gabadadze, and A. Vainshtein, Phys. Rev. D 65, 044026 (2002), eprint arXiv:hepth/0106001.
 H. Oyaizu, Phys. Rev. D 78, 123523 (2008), eprint 0807.2449.
 R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (Computer Simulation Using Particles, New York: McGrawHill, 1981, 1981).
 A. Klypin and J. Holtzman, ArXiv Astrophysics eprints (1997), eprint astroph/9712217.
 A. Brandt, Math. Comput. 31, 333 (1977).
 A. Brandt, Tech. Rep. (GMDStudien) Nr. 85, Gesellschaft für Mathematik und Datenverarbeitung, Schloss Birlinghoven, Sankt Augustin, Germany (1984).
 W. L. Briggs, V. E. Henson, and S. F. McCormick, A multigrid tutorial (2nd ed.) (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2000), ISBN 0898714621.
 W. Fang, S. Wang, W. Hu, Z. Haiman, L. Hui, and M. May, Phys. Rev. D 78, 103509 (2008c), eprint 0808.2208.
 A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), eprint astroph/9911177.
 R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, MNRAS 341, 1311 (2003), eprint arXiv:astroph/0207664.
 J. L. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. S. Warren, G. Yepes, S. Gottlober, and D. E. Holz, ArXiv eprints 803 (2008), eprint 0803.2706.
 I. Laszlo and R. Bean, Phys. Rev. D 77, 024048 (2008), eprint 0709.0307.
 S. D. M. White, J. F. Navarro, A. E. Evrard, and C. S. Frenk, Nature (London) 366, 429 (1993).
 V. R. Eke, S. Cole, C. S. Frenk, and J. Patrick Henry, MNRAS 298, 1145 (1998), eprint arXiv:astroph/9802350.
 S. Borgani, P. Rosati, P. Tozzi, S. A. Stanford, P. R. Eisenhardt, C. Lidman, B. Holden, R. Della Ceca, C. Norman, and G. Squires, Astrophys. J. 561, 13 (2001), eprint arXiv:astroph/0106428.
 A. Vikhlinin, A. V. Kravtsov, R. A. Burenin, H. Ebeling, W. R. Forman, A. Hornstrup, C. Jones, S. S. Murray, D. Nagai, H. Quintana, et al., Astrophys. J. 692, 1060 (2009), eprint 0812.2720.
 A. Lue, R. Scoccimarro, and G. D. Starkman, Phys. Rev. D 69, 124015 (2004), eprint arXiv:astroph/0401515.
 D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. Blanton, R. C. Nichol, R. Scranton, H.J. Seo, M. Tegmark, Z. Zheng, et al., Astrophys. J. 633, 560 (2005), eprint arXiv:astroph/0501171.
 W. J. Percival, S. Cole, D. J. Eisenstein, R. C. Nichol, J. A. Peacock, A. C. Pope, and A. S. Szalay, MNRAS 381, 1053 (2007), eprint 0705.3323.
 W. Hu, private communication.
 L. Fu, E. Semboloni, H. Hoekstra, M. Kilbinger, L. van Waerbeke, I. Tereno, Y. Mellier, C. Heymans, J. Coupon, K. Benabed, et al., A&A 479, 9 (2008), eprint 0712.0884.
 URL http://fftw.org/.
 S. F. Shandarin, Astrophysics 16, 439 (1980).
 A. V. Kravtsov, A. A. Klypin, and A. M. Khokhlov, Astrophys. J. Supplement 111, 73 (1997), eprint arXiv:astroph/9701195.
 W. Hu and I. Sawicki, Phys. Rev. D 76, 104043 (2007), eprint arXiv:0708.1190.