# Efficient simulations of large scale structure in modified gravity cosmologies with comoving Lagrangian acceleration

###### Abstract

We implement an adaptation of the COLA approach, a hybrid scheme that combines Lagrangian perturbation theory with an N-body approach, to model non-linear collapse in chameleon and symmetron modified gravity models. Gravitational screening is modeled effectively through the attachment of a suppression factor to the linearized Klein-Gordon equations.

The adapted COLA approach is benchmarked, with respect to an N-body code both for the CDM scenario and for the modified gravity theories. It is found to perform well in the estimation of the dark matter power spectra, with consistency of 1% to h/Mpc. Redshift space distortions are shown to be effectively modeled through a Lorentzian parameterization with a velocity dispersion fit to the data. We find that COLA performs less well in predicting the halo mass functions, but has consistency, within uncertainties of our simulations, in the relative changes to the mass function induced by the modified gravity models relative to CDM.

The results demonstrate that COLA, proposed to enable accurate and efficient, non-linear predictions for CDM, can be effectively applied to a wider set of cosmological scenarios, with intriguing properties, for which clustering behavior needs to be understood for upcoming surveys such as LSST, DESI, Euclid and WFIRST.

## I Introduction

The nature of the unknown mechanism responsible for the accelerated expansion of the universe, as measured by Type 1a supernovae Perlmutter et al. (1999); Riess et al. (2004), baryon acoustic oscillations (BAO) in galaxy clustering Eisenstein et al. (2005); Percival et al. (2007, 2009); Kazin et al. (2014), and the Cosmic Microwave Background (CMB) Spergel et al. (2003); Ade et al. (2013, 2016), commonly labeled as “Dark Energy”, is one of the most challenging, open questions in modern cosmology. Assuming that Einstein’s General Relativity (GR) is the correct framework to describe gravity at large scales, the recent accelerative phase can be driven either by a cosmological constant term with negative pressure, or by introducing a scalar field called quintessence Wetterich (1988); Ratra and Peebles (1988); Copeland et al. (2006). The necessary value of to account for the observed acceleration rate is extremely small, however, when compared to the predictions from high energy physics, and as a result it has to be fine-tuned Weinberg (1989). Such an unattractive feature, together with the need to explore all other alternatives, has motivated the development of theories in which GR breaks down at large scales Carroll et al. (2004); Capozziello et al. (2003); Nojiri and Odintsov (2006), the so called modified gravity (MG) theories. The existence of MG theories with massive scalar fields can reproduce the recent accelerative phase, however, the fifth forces that arise as a result of their coupling to matter would, in principle, cause large deviations from the tight experimental constraints of GR in the solar system Will (2006).

As a consequence, MG models can only be viable if they reduce to the successful GR phenomenology in the local dense environments (eg. Earth, Solar System) through a restoring screening mechanism Khoury (2010, 2013). Based on the qualitative features of the screening mechanism they exhibit, such schemes are commonly classified in various broad classes: the “chameleons” Khoury and Weltman (2004a, b), where the scalar fields become massive and decouple in regions of high Newtonian potential, the kinetic/“-Mouflage” models Babichev et al. (2009); Dvali et al. (2011), in which the deviations are screened when fifth forces exceed some critical model dependent value and the Vainshtein mechanism Vainshtein (1972), that reproduces GR when large derivatives of the fifth forces are experienced. Similar, in terms of phenomenology, with the chameleons are the symmetrons Hinterbichler and Khoury (2010); Olive and Pospelov (2008), which exhibit the additional property of a vanishing coupling in dense regions through symmetry restoration.

This rich spectrum of MG models are theoretically viable and offer observational consequences that are potentially distinguishable from CDM through a variety of astrophysical characteristics.
A number of spectroscopic and photometric Large Scale Structure (LSS) surveys both currently underway, e.g. Dark Energy Survey (DES) (Abbott et al., 2005), HyperSuprimeCam (HSC) ^{1}^{1}1http://sumire.ipmu.jp/en/, and eBOSS (Comparat et al., 2012), and coming online in the coming decade, e.g. DESI Levi et al. (2013), PSF, LSST Abell et al. (2009), Euclid Laureijs et al. (2011) and WFIRST Spergel et al. (2013), will probe the properties of gravity with remarkable precision in both the linear and non-linear regimes, using galaxy clustering, cluster counts, gravitational lensing and peculiar velocities. They offer an unprecedented opportunity to test the landscape of modified gravity theories observationally with respect to the simplest, CDM scenario. As a result, simulating the structure formation in the the linear, mildly non-linear and non-linear regimes is necessary for both CDM and all the alternatives.

A variety of analytical, semi-analytical and numerical approaches have been used to study CDM and dark energy scenarios in the non-linear regime. Lagrangian perturbative techniques up to first Zeldovich (1970); Peebles (1980) or second order Bouchet et al. (1995), have been shown to produce accurate results for CDM in the linear and mildly non-linear scales without having to perform a complete numerical treatment of structure formation. They fail to achieve the desired accuracy, however, at smaller, non-linear scales for which a full N-body simulation is required. In light of the computational resources necessary for N-body simulations, and given the successes of Lagrangian Perturbation Theory (LPT), hybrid schemes have been proposed, with the aim of combining the strengths of both approaches. In this paper we focus on the Comoving Lagrangian Acceleration (COLA) hybridization scheme Tassev et al. (2013). By evolving the large scales analytically using LPT and the small scales exactly with a full N-body treatment, the COLA method manages to produce accurate results deep in the mildly non-linear regime with only a few number of time steps, making it possible to produce fast results in exchange for some accuracy.

In modified gravity simulations, the need to accurately capture the effects of the fifth forces and the screening mechanism adds a new layer of complexity. For an exact description, one needs to solve the full Klein-Gordon equation, whose non-linearities render the procedure both challenging and computationally expensive. It is natural consequently to investigate whether an inexpensive, approximate scheme can be used instead. A linear treatment of the perturbation equations, together with the linearized Klein-Gordon equation it produces may seem efficient at first, but a more careful examination shows Brax et al. (2012, 2013) that it fails to incorporate the non-linear screening effects and gives poor results. Effective approaches Winther and Ferreira (2015) have managed to implement screening successfully, however, following a phenomenological path. An ineffective but computationally fast linearized scheme, can be combined with the attachment of a screening factor for a spherically symmetric configuration, to speed up MG simulations without the sacrifice of much accuracy.

Given the success of Lagrangian approaches in CDM simulations and the need to develop effecient, but representative, realizations of the LSS in different cosmological scenarios, it is natural to see alternative routes in MG models. The benefits of LPT have already been discussed in the context of generating initial conditions, appropriate for coupled scalar field cosmologies Li and Barrow (2011) or MG models Valkenburg and Hu (2015) . In this paper, we study the effectiveness of the COLA hybrid scheme, in which the linear scales are evolved exactly using LPT and the non-linear ones using N-body simulations, for MG scenarios. As far as the N-body component is concerned, the fifth force calculation lies in the solution of the linearized KG equation and an approximate screening implementation through the thin shell factor for a dense sphere, similar to Winther and Ferreira (2015). In chameleon-type (and symmetron) models, a scalar field acquires a very large mass within a massive object and consequently decouples due to the Yukawa suppression, so essentially only a fraction of the total mass (thin shell) contributes to the fifth force.

The layout of the paper is as follows: in Sec. II we first review the MG models studied and the non-linear approaches used in the analysis. In Sec. III we present our results, assessing the performance for the scheme to predict a number of LSS observables, including the matter power spectrum, the redshift space distortions, and halo mass function, before summarizing the findings and discussing implications for future work in Sec. IV.

## Ii Formalism

### ii.1 Modified gravity and screening models

A wide class of viable scalar-tensor theories have been shown to be described by a Horndeski Lagrangian Horndeski (1974); Deffayet et al. (2011). Using a general single scalar field Lagrangian, in the Einstein frame, written in terms of a scalar field and its derivatives,

(1) |

where R is the Ricci scalar, the scalar field, the reduced planck mass and is the Lagrangian for the matter sector, in which the matter fields are non-minimally coupled to the scalar field with a dimensionless coupling constant . In the chameleon and symmetron models, the properties of the single scalar field can be described by a simple, scalar field Lagrangian

(2) |

where is the self-interacting potential. Varying the action gives us the equations of motion for the scalar field, the Klein-Gordon equation

(3) |

where the effective potential combining the self-interaction potential and coupling term is given by

(4) |

The chameleon screening mechanism lies in the fact that the effective mass of the scalar field calculated at the minimum, , which is given by

(5) |

has to be positive. For the chameleon theories, this requirement is guaranteed through the interplay between a monotonically decreasing potential and an increasing coupling. In the symmetron model, on the other hand, the viability is restored using a “Mexican hat” symmetry breaking potential Hinterbichler and Khoury (2010), the behavior of which still gives rise to a positive density-dependent mass.

The observational consequences of such models can be demonstrated by extracting the scalar field profile, , produced by the density profile

(6) |

where is the radial distance from the center of a compact spherically symmetric configuration of density and radius (not to be confused with the Ricci scalar R), that is isolated on a uniform density background . Under spherical symmetry, (3) becomes

(7) |

Even though (7) does not have, in principle, an analytical solution, accurate approximations can be performed for two different configurations, that correspond to opposite regimes with respect to screening Khoury and Weltman (2004a, b). The first case is that of a large, strongly perturbing object of very large density , for which the interior field is forced to acquire the value that corresponds to the minimum of the effective potential, and the scalar field profile outside the object is given by

(8) |

The corresponding fifth-force experienced by a unit mass particle outside the object is

(9) |

where are respectively the background values of the mass and coupling and M the mass of the object. Given the characteristic large values of the Compton wavelength , the scalar field is essentially free within our scales of interest and the Yukawa suppression can be neglected in (9),

(10) |

The above approximation is valid when the “screening factor” is

(11) |

which also defines the criterion for the existence of a thin shell Hinterbichler and Khoury (2010); Khoury (2013), whose mass is the fraction of the total that actually contributes to the fifth force, due to the strong Yukawa suppression deep inside dense objects. The Newtonian gravitational potential is denoted by in (11). On the other hand, when linear perturbation theory is valid, which is the case when , the linearized form of (7) gives

(12) |

for the fifth force. Based on (10)-(12), we see that in the linear regime the fifth force is the same as the Newtonian force with a coupling and deep in the non-linear (screened) regime, it is suppressed by the thin shell factor (11).

Furthermore, it should be also noted that, as shown in Brax et al. (2012), one can derive a pair of functions , for the characterization of a model within the above framework. Unlike models with constant couplings, symmetrons exhibit an additional form of screening Hinterbichler and Khoury (2010); Olive and Pospelov (2008) due to the fact that in dense environments symmetry is restored and the coupling vanishes.

Adopting this formulation, linear perturbation theory gives Li and Zhao (2009); Brax et al. (2012) for the growth of CDM density perturbations in the quasi-static limit and for sub-horizon scales

(13) |

with

(14) |

where is scale factor, with today, and is the comoving wavenumber.

The effects of gravity modifications at the linear approximation are incorporated in the second term. For very large scales and/or early times (GR regime), and (13) reduces to the standard GR expression in the weak gravity regime, where the Newtonian gravitational potential is given by the Poisson equation,

(15) |

When however (scalar-tensor regime), the second term becomes significant and gives the linearized Klein-Gordon equation for the fifth potential

(16) |

with the real space expression being

(17) |

#### ii.1.1 The model

theories Carroll et al. (2004) are widely-studied modified gravity scenarios, that give rise to acceleration on cosmic scales and can be incorporated Brax et al. (2008) into the chameleon formalism with a constant coupling . The first model we tested thus, was the Hu-Sawicky model Hu and Sawicki (2007) with a scalar field mass

(18) |

where

(19) |

where is the Hubble Constant and and are, respectively, the dark energy and dark matter fractional energy densities today. The mass takes the form

(20) |

Furthermore, the screening factor is given by

(21) |

and are the model’s free parameters. In this paper, we consider the model for and . These describe cosmologically viable scenarios whose non-linear properties have been simulated using the full Klein-Gordon equation Zhao et al. (2011); Winther and Ferreira (2015) with which our results can be compared.

#### ii.1.2 The symmetron model

The general framework laid previously, can also incorporate the symmetron model, with a “Mexican hat” symmetry breaking potential Hinterbichler and Khoury (2010), for which scalar fields couple to matter after , with

(22) |

and the coupling vanishes for , when symmetry is restored. The screening factor for this model becomes (Davis et al., 2012; Winther and Ferreira, 2015)

(23) |

We consider this model with values and which again have been shown Davis et al. (2012) to predict deviations consistent with experimental constraints. It should be also pointed out that, as explained previously, models of this type exhibit field dependent couplings which cause additional screening due to the coupling suppression in dense environments, where symmetry is again restored. This effect is not taken into account in our approximate scheme.

### ii.2 Simulating non-linear clustering

#### ii.2.1 The N-body method

The COLA code has been loosely based on A. Klypin’s PM code Klypin and Holtzman (1997), and this motivates the latter’s use as a comparison for our approximate scheme’s effectiveness. It is also a simple and representative implementation of a Particle-Mesh (PM) N-body code. N-body simulations for MG using the PM code have been performed previously Stabenau and Jain (2006); Laszlo and Bean (2008); Khoury and Wyman (2009). For each scenario, we consider 10 simulated realizations, initialized at an initial redshift , at which density perturbations on the scales we study are linear. After providing a linear power spectrum from the cosmological code CAMB Lewis et al. (2000) for the desired CDM cosmology at the time , particles are placed in our simulation box with side L=200 Mpc/h, in a mesh of , using 1st order Lagrangian Perturbation Theory (Zel’dovich approximation) Zeldovich (1970). The parameters that define our background CDM cosmology are , , , and . The particle positions are updated, using 500 time steps, through the displacement equation:

(24) |

In Fig. 1, it is shown that the choice of 500 iterations, which corresponds to steps of in the scale factor, guarantees convergence at the 0.08% level.

In MG cosmologies, the modified geodesic equation gives, in the weak gravity regime, the modified version of (24),

(25) |

where the term is negligible given observational constraints from variations of constants Winther and Ferreira (2015). Equation (25), which also holds for the full non-linear KG description, forms a closed system of equations with (15) and (17) that are solved in the Fourier space for the potentials and .

The linearized form of KG equation, (17), does not incorporate the screening effects. To account for the screening effect, we adopt an effective parameterization similar to the one proposed in Winther and Ferreira (2015). In section Sec. II.1, we showed that the linear solution for the fifth force, (10), is suppressed by the screening factor deep in the non-linear regime. As a result, we incorporate the screening effects by explicitly attaching the screening factor to the fifth force in accordance with (10)-(12) and (25),

(26) |

To interpolate properly between the screened and the unscreened regime we set

(27) |

Within our approximate scheme, the functions and have been set equal to the background ones and correspondingly, which has been shown to be a good approximation in (Winther and Ferreira, 2015).

#### ii.2.2 The COLA method

The fact that N-body codes manage to simulate the Large Scale Structure accurately but at a significant computational cost, has motivated the development of several analytical perturbative techniques to avoid a full blown N-body simulation. Lagrangian Perturbation Theory (LPT) Zeldovich (1970); Bouchet et al. (1995) works perturbatively in a Lagrangian displacement field and manages to give accurate results in the Linear and the Mildly Non-Linear regime. However, it quickly fails to capture the non-linearities associated with the smaller scales and consequently it underestimates significantly the power at large . Given that we have to choose between accurate, but expensive N-body simulations and fast but approximate perturbative techniques, it is reasonable to ask whether one can efficiently combine the benefits of both approaches. Such a hybrid method, named COmoving Lagrangian Acceleration (COLA) was proposed in Tassev et al. (2013). Here we outline the basic framework and its modifications for MG, while details can be found at Tassev et al. (2013). The particle comoving positions are decomposed as a sum of two pieces, in the “manifestly” exact form

(28) |

By defining a new time variable , where is conformal time, (24) can be cast in the simpler form

(29) |

with , and . In the Lagrangian description , with the initial Eulerian position and the Lagrangian displacement and

(30) |

One can now solve for the residual displacement in CDM

(31) |

where and are the first and second order growth factors, respectively, and , are the Zel’dovich and second order LPT displacements. The fact that the LPT piece is evolved analytically and we only solve numerically for , can be interpreted as working on a frame that is co-moving with observers that follow LPT trajectories.

can be discretized using a Leapfrog scheme Quinn et al. (1997) to get the core COLA equations for each particle’s position and velocity change between the times

(32) |

where a tilde denotes a quantity in units of .

Initial conditions are produced using the 2LPT initial conditions code (2LPTic) Scoccimarro (1998) which does so by performing LPT up to second order. In a CDM cosmology, growth functions and are scale independent Zeldovich (1970); Bouchet et al. (1995) and one only needs to produce an LPT snapshot for for both generating initial conditions and obtaining the LPT terms at the different timesteps. In such a case, the LPT displacements are given by , and (32) reduces to the standard COLA CDM scheme (with the fifth force term omitted). Initial conditions and background cosmology are produced, for 10 realizations, for the same cosmological parameters as used in the PM code, at the initial redshift z=9.0 which has been shown Tassev et al. (2013) to work well for COLA in CDM. The simulation box size, number of particles and mesh size are the same as used in the PM code. It should be noted though that we don’t perform a comparison of the codes by initiating both with identically seeded initial conditions, but instead, we compare the statistical consistency of the means of the 10 runs for each of the two techniques with the respective sets each using different random generated seeds. In its initial formulation, COLA was used with 10 time steps, which enables accurate predictions down to h/Mpc, which can be also seen in Fig. 1, where the CDM power spectrum by COLA is presented for various choices of time steps. By increasing the number of steps to 50, still significantly fewer than the typical number of iterations performed in a standard N-body code, we can provide accuracy down to smaller scales, h/Mpc. The CDM COLA run-time in this set up is 10 times shorter than that of the PM code. COLA’s accuracy as a function of the number of time steps used, is further discussed in section III.1.

When gravity is modified, the core equations need to be changed appropriately to account for the additional fifth forces and the screening effects. As in the N-body code, one can use an approximate framework to model the modified gravity effects both on the growth rate and screening: solving the linearized KG equation (16) and attaching the screening factor (27) to fifth force term in (32). For the LPT component of COLA, one must consider that, in MG theories, the growth factors and become scale dependent.

In Figure 2, we summarize the first and second order growth factors for a chameleon and symmetron scenario. In each case as we approach late times, for the chameleon model, and for the symmetron, we find significant scale dependent deviations from CDM at the level of 10% for and 25% for at k=0.1 h/Mpc today, which means that not all Fourier modes evolve the same way with time Valkenburg and Hu (2015). This causes the LPT trajectories of a given particle to bend, in principle. As a consequence one has to be very cautious about how to obtain the LPT terms at the different times. We briefly outline the application of LPT to scale dependent growth functions in MG in appendix A.

Unlike the CDM case, here the growth factors’ scale dependent nature does not allow one to evolve the Zel’dovich and order displacements with a single scale independent function for all scales. To account for that, we have considered two alternative modifications to COLA. In the first approach, we create an MG version of COLA that calculates the LPT displacements numerically at each time step using an MG version of 2LPTic. The relevant LPT terms in (32) are calculated after Fourier transforming (56) and (57). Besides the modified N-body component, the fact that we have to solve numerically for the Lagrangian terms at every discrete time step increases the computational cost significantly. In a second approach, we utilize the fact that the LPT part of the scheme serves to evolve the linear scales, for which the MG deviations with respect to CDM are known to be small for most times and adopt an approximate scheme in which only the N-body part is modified and the CDM solutions are used for the Lagrangian displacements. The resulting scheme has the same N-body component as in (32) and the known CDM LPT terms, in which the Lagrangian displacements are evolved with and .

(33) |

A comparison of the two approaches for the model with , for which we expect the largest modifications, is shown in Fig. 3. One comparison tracks a given particle inside our volume during the simulation and we also compare the resulting power spectra using both schemes. We find excellent agreement between the approximate and fully modified schemes, in both the linear and mildly non-linear regimes. We also find very small differences for the position and displacement vectors (magnitude & direction) with differences in angular orientation of at most arcseconds, and differences in the magnitude of steps less than , which result in power spectra that have a fractional difference no larger than 0.3% today. The approximate scheme takes under half the run time of the full implementation. In light of these results, we adopt the approximate scheme in the COLA simulations used in this analysis. This has the great advantage of not having to solve numerically for the LPT displacements at every time step, without sacrificing much accuracy.

## Iii Analysis/Results

### iii.1 Modified gravity results

In this section we present the results of the assessment of COLA’s performance with respect to the predicted power spectra, redshift space distortions (RSD) and dark matter halos for the modified gravity scenarios and CDM. For every given model and choice of parameters, simulations have been performed using both COLA and the PM code.

#### iii.1.1 Power spectra

To appropriately benchmark the COLA performance for modified gravity, we first compare the performance of COLA for CDM. In Fig. 4, we show the CDM power spectra as obtained by both codes, together in comparison with the fit by Smith et al. (2003). The two results agree well within a standard deviation of each other for all scales, but start to, underestimate power, consistently with one another, by h/Mpc, relative to higher resolution simulations. For that reason, we choose to compare performance down to a scales with h/Mpc, while the Nyquist wavenumber, for our simulation, is h/Mpc.

In Fig. 5, the fractional difference in the power spectra is plotted for all our models and both codes are found to agree with each other well within one standard deviation, with the differences being smaller than 1%. Our results demonstrate the consistency between COLA and the N-body approach using the approximate scheme. In turn this connects with previous work that has shown, in general, the good degree of consistency of this approximate scheme with N-body simulations using the full Klein Gordon for the same models in the literature Zhao et al. (2011); Davis et al. (2012); Brax et al. (2012); Winther and Ferreira (2015). In particular the results for the & models are in excellent agreement with the literature for all scales. Our results confirm findings in Winther and Ferreira (2015), in studying the effectiveness of the linearized screening schema: for the lowest screening model, with , and the symmetron model the effective screening parameterization, respectively, under and overestimates the power, relative to the full KG simulation, at the non-linear scales.

Fig. 6 shows our COLA scheme’s accuracy in predicting the fractional difference in the power spectra for the highest deviation model, , as a function of the number of time steps used, for one realization. We find that using 50 time steps provides excellent convergence, at the level of 0.9%, to the scales we want to consider, Mpc/h. Using 30 time steps provides convergence at the level of 8% at Mpc/h.

#### iii.1.2 Redshift space distortions

A great amount of observational effort is being invested in studying the three-dimensional Large Scale Structure (LSS) through spectroscopic galaxy surveys that measure precise redshifts. Among various challenges faced by such measurements, the observed clustering structures appear distorted in redshift space.

Density perturbations give rise to peculiar velocities with respect to the Hubble flow, which result in the redshift space position , being different than the real space position , with the relationship between them taking the form

(34) |

By we denote the peculiar velocity and by the unit vector along the line of sight. At linear scales, coherent motions of galaxies that tend to collapse within an overdense region, cause it to appear squashed in redshift space. As shown by Kaiser (1987), in the distant observer approximation, such an overdensity will be distorted in the redshift space:

(35) |

where is the angle between the peculiar velocity and the line of sight in space, , and the linear growth rate,

(36) |

with subscript ‘’ to differentiate it from the function. Such an effect gives rise to, based on (35), an overestimation of the power spectrum measured in the redshift space:

(37) |

where we introduced the factor (not to be confused with the coupling ) to account for the galaxy bias , with for cold dark matter. Averaging (37) over all directions, gives the order piece

(38) |

At smaller, non-linear, scales the random incoherent velocities of galaxies within virialized structures cause overdense regions to appear elongated along the line of sight (“Fingers of God”), causing suppression of power. An exact quantitative treatment of the phenomenon is hard, due to the complicated nature of the small-scale velocity correlations and as a result phenomenological approaches have been proposed. Such models Peacock and West (1992) treat the line-of-sight distortion as a radial convolution of the correlation function (including the Kaiser boost) with an incoherent velocity distribution

(39) |

where and are the perpendicular and parallel components. Assuming a Gaussian velocity distribution Peacock and West (1992), the Fourier space expression would then be

(40) |

with being the comoving distance dispersion that is related Wang et al. (2013) to the velocity dispersion through

(41) |

Even though the exponential term in (40) is reasonable as a damping term for capturing the non-linear power suppressions, it has been noted Davis and Peebles (1983) that an exponential pairwise velocity distribution

(42) |

is a better fit. This gives rise to the dispersion model Peacock and Dodds (1994)

(43) |

in which the damping effects are incorporated through a Lorentzian term and (or ) is considered a free parameter to be fitted to the data. It should be noted that is actually scale and bias dependent, which is one of the limitations the dispersion model faces Scoccimarro (2004). The above description can still prove to be a very useful tool for obtaining an effective non-linear velocity dispersion parameter and thus quantifying the non-linear FoG effect. Integrating (43) over all directions gives the monopole piece, which can be fitted over the results to obtain . This is slightly different than other approaches: Angulo et al. (2008) proposed attaching a simple factor to the Kaiser boost with being a free parameter, loosely related to , while Kwan et al. (2012) suggested attaching a free function and marginalized over the parameters , and for to account for the uncertainties in constraining the effects of modified gravity on the RSD power spectrum.

Analytic prediction | ||||
---|---|---|---|---|

Scenario | ||||

CDM | 1.35 | 0.462 | 0.556 | 567 |

, | 1.35 | 0.463 | 0.555 | 605 |

, | 1.38 | 0.491 | 0.512 | 714 |

, | 1.42 | 0.541 | 0.443 | 834 |

Symmetron | 1.35 | 0.464 | 0.554 | 611 |

Through the mapping (34), we obtained redshift space power spectra for all of the simulated models, with the results presented in Fig.7. We compared the large scale results to analytic predictions arising from the linear growth rate, and also used the monopole model in (43) to obtain an effective velocity dispersion damping factor for the FoG effect. The results are summarized in Fig. 7 and Table 1.

We first benchmarked COLA’s performance for CDM. We see that the PM and COLA codes’ RSD predictions for CDM do not differ by more than 0.5% at all the scales of interest and agree remarkably well with the analytical prediction, with expected values of =0.467 and =1.354, assuming , with . At smaller scales , the “Fingers of God” effect quickly dominates, and causes power suppression and find this suppression is well modeled by (43) with .

For the MG models, the additional fifth forces cause the redshift space power spectra to have, in principle, different shapes. In large scales, the enhanced clustering results in higher coherent velocities of collapse into overdense regions which translates to a higher boost in the RSD power spectra with respect to GR, translating into higher values for the growth rate and a lower . For lower magnitude modifications, the suppression of the fifth forces gives results that tend to the CDM prediction. At smaller scales, the fifth forces cause higher random velocity dispersions inside virialized structures, making the damping effects stronger in MG.

These combined effects cause the redshift space distortions to be more pronounced in MG compared to GR. This can be clearly seen in the upper left panel in Fig.(7). As expected, the redshift space distortions vary from the most pronounced, in the lowest screening model, to very small deviations from GR in the strong screening regime. For the same models, the redshift space power spectra from the PM code agree with COLA well within a standard deviation. The results using the approximate schema are in good agreement with full non-linear MG N-body simulations for redshift space distortions in gravity performed by Jennings et al. (2012). For all the models, COLA predicts deviations that are 0.5% more pronounced (higher in large scales, smaller in small scales) than the PM code.

#### iii.1.3 Halo Mass Function

To determine the halo mass function we identify halos in the simulations using the Rockstar halo finder Behroozi et al. (2013) for all models. In Fig. 8 we show the comparison of the halo mass function predicted by COLA and the PM code, together with a high accuracy result by Murray et al. (2013). COLA and PM are found to be in a better than 2.5% agreement in the lower and intermediate mass range, while in the highest mass bins there is a maximum difference of 10%.

In Fig. 9, we plot the fractional difference in the halo mass function with respect to CDM, for all of our models. The COLA and PM code results agree in general, well within the standard deviation from the averaged suite of simulations. In particular, In the and models, the PM code predicts a fractional boost in the halo mass function that is higher than COLA’s by 2% and 2.5%, for the lower and intermediate bins, while in the highest bin COLA gives a boost larger by 5% and 3% correspondingly. For the and symmetron models, the differences are 1% and smaller, with the PM code giving greater number counts for the two mass bins below and COLA being higher for the bin over . The differences between the predictions in each case and especially in the high mass bin, are within, and likely largely resulting from, the differences observed in the CDM benchmarking of the mass functions.

While we do not perform a simulation with the full non-linear Klein-Gordon equation, we note that compared to other full KG treatments in the literature Winther and Ferreira (2015); Zhao et al. (2011), our method performs well and only slightly underestimates the mass function for the & models, in accordance with the general features noticed in the power spectra discussion. In agreement with Winther and Ferreira (2015), we observe an underestimation of halos in the lower end of our mass range (around M ) for the model, indicating too much screening, and an overestimation of the mass function for the symmetron model.

## Iv Conclusions

In this paper, we have implemented a hybrid scheme, that combines Lagrangian Perturbation Theory and N-body approaches, to numerically characterize the evolution of large scale structure in chameleon and symmetron, modified gravity, theories which exhibit gravitational screening in the non-linear regime. LPT is used to evolve linear scales analytically in combination with a full N-body approach that is used for the non-linear scales to reduce computational costs. An effective screening scheme is implemented in place of a solution to the full Klein-Gordon equation for the fifth potential, in which an effective suppression factor is attached to the real-space linearized perturbations.

We demonstrate that while in MG spatial modes evolve differently in LPT (and can have deviations from the nominal GR geodesic paths), the scheme can be further simplified, for the models we studied, by using a displacement coordinate system based on scale-independent CDM growing modes combined with a modified, screened Poisson equation. We note that, while this approximate scheme works well for the chameleon and symmetron models we consider, it should always be tested against the exact LPT solution for a new modified gravity model.

Our method was applied on the and symmetron models and it was tested against power spectra, redshift space distortions and dark matter halo mass functions, using a fiducial number of 50 time steps. At the same time, we assessed our hybrid’s performance against simulations from a pure N-body code with the same screening implementation for the same models, using 500 iterations.

With regards to power spectra, we found COLA to be in better than 1% agreement with the N-body code at all scales for all the models studied. Note that the effective screening scheme we use has previously been shown to be in good agreement with results using the full non-linear Klein Gordon in an N-body implementation Winther and Ferreira (2015). We find, as was discussed in (Winther and Ferreira, 2015), that the effective screening approach does underestimate power, relative to that found in solving the full Klein-Gordon Zhao et al. (2011), as one moves into the fully non-linear regime (), however this is also beyond the regime of applicability of COLA’s scheme.

COLA and the N-body code are in better than 0.5% agreement with respect to redshift space distortions for all the scales and models of interest. The distortions were modeled by attaching the linear Kaiser factor for the enhancement at large scales and a Lorentzian dispersion factor for the small scale suppression due to incoherent motions within virialized structures. We find that the monopole is a well fit using an effective pairwise velocity dispersion as a fitting parameter to quantify the suppressions at non-linear scales. The additional fifth forces present in the chameleon and symmetron models, cause the redshift space distortions to be more pronounced with respect to CDM. This can be seen by the larger boosts in linear scales due to the higher coherent velocities, and by the stronger suppressions in the non-linear scales because of the higher values of the velocity dispersion. The adapted COLA scheme gives reasonable results for the predicted fractional boost in the halo mass function relative to CDM, with the differences between the N-body and COLA results in the halo mass function estimation most likely being due to the difference between the two codes in CDM.

In this paper, we have focused on chameleon and symmetron-type scalar-tensor theories, but it would be very interesting to see how well this scheme performs for the simulation of other screening mechanisms as well such as the Vainshtein mechanism (Vainshtein, 1972), as well as other dark energy models, such as those with non-minimal couplings between dark matter and a quintessence scalar field (Li and Barrow, 2011). Given the level of consistency between COLA and the N-body predictions for the monopole of the redshift power spectrum, it would also be interesting to investigate the COLA scheme’s ability to capture higher order moments of the angular power spectrum to, for example, calculate the ratio of the quadrupole to monopole moments to estimate in a way that is robust to systematic effects from incomplete modeling of the nonlinear distortions Landy and Szalay (2002); Jennings et al. (2012).

Many theories being considered as explanations for cosmic acceleration have tantalizing predictions in the non-linear regime but also present computational challenges in modeling them. With a suite of next-generation large scale structure surveys, including LSST, DESI, Euclid and WFIRST, starting in next few years, there is an unprecedented opportunity to measure the properties of large scale structure clustering as it transitions from linear to mildly and then strongly non-linear scales, and using multiple tracers. The results presented here demonstrate that COLA, proposed to enable accurate and efficient, non-linear predictions for CDM, is a viable approach to study non-linear collapse for a broader portfolio of cosmological scenarios. For example, in work that has followed our paper in Winther et al. (2017), the effectiveness of the COLA approach has also been studied in the and nDGP models, and was shown to perform very well in predicting the fractional deviations with respect to the CDM power spectra and halo mass functions, using a small number of time steps.

## Acknowledgments

We want to thank Hans Winther and Pedro Ferreira for useful comments on our paper and Risa Wechsler for helpful discussions on halo mass function performance with the COLA algorithm. We would also like to thank an anonymous referee, whose valuable comments helped improve and clarify this manuscript. The work of GV and RB is supported by NASA ATP grant NNX14AH53G, NASA ROSES grant 12-EUCLID12- 0004 and DoE grant DE-DE-SC0011838.

##
Appendix A

Lagrangian perturbation theory in Modified Gravity

LPT Bouchet et al. (1995) works perturbatively in a displacement field

(44) |

where and are the initial and final comoving Eulerian particle positions. In this formulation, all the information is reflected in the mapping through the displacement field. Working up to first order gives the so-called Zel’dovich approximation in CDM, for which the can be decomposed into a product of temporal and spatial factors

(45) |

and

(46) |

with being the Gaussian density field generated by an initial linear power spectrum and the scale independent first order growth factor, given by

(47) |

In an MG scenario, the growth factor is not scale independent any more and

(48) |

where

(49) |

in Fourier space. This implies that particle trajectories, unlike in CDM, are not straight lines Valkenburg and Hu (2015). (46) and (48) indeed give

(50) |

The second term, responsible for the trajectory bending, vanishes when the growing mode is scale independent, in which case one recovers the standard CDM Zel’dovich approximation. When working up to second order (2LPT), we have in a similar fashion

(51) |

where the second order growth factor is given by

(52) |

For the early times, the spatial part is given by

(53) |

In the MG case, we will have again

(54) |

with the scale dependent second order growth factor that obeys

(55) |

in the Fourier space. The fact that all of our models recover GR at early times, guarantees that the early time spatial part is still given by (53). In our implementation of the full MG COLA scheme, a suitably modified version of 2LPTic produces the LPT terms at every time step, through the Fourier space versions of (48) and (54)

(56) |

and also the same for the accelerations

(57) |

## References

- Perlmutter et al. (1999) S. Perlmutter et al. (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), eprint astro-ph/9812133.
- Riess et al. (2004) A. G. Riess et al. (Supernova Search Team), Astrophys. J. 607, 665 (2004), eprint astro-ph/0402512.
- Eisenstein et al. (2005) D. J. Eisenstein et al. (SDSS Collaboration), Astrophys.J. 633, 560 (2005), eprint astro-ph/0501171.
- Percival et al. (2007) W. J. Percival et al., Mon. Not. Roy. Astron. Soc. 381, 1053 (2007), eprint 0705.3323.
- Percival et al. (2009) W. J. Percival et al. (2009), eprint 0907.1660.
- Kazin et al. (2014) E. A. Kazin, J. Koda, C. Blake, and N. Padmanabhan (2014), eprint 1401.0358.
- Spergel et al. (2003) D. N. Spergel et al. (WMAP), Astrophys. J. Suppl. 148, 175 (2003), eprint astro-ph/0302209.
- Ade et al. (2013) P. Ade et al. (Planck Collaboration) (2013), eprint 1303.5076.
- Ade et al. (2016) P. A. R. Ade et al. (Planck), Astron. Astrophys. 594, A13 (2016), eprint 1502.01589.
- Wetterich (1988) C. Wetterich, Nuclear Physics B 302, 668 (1988).
- Ratra and Peebles (1988) B. Ratra and P. J. E. Peebles, Phys. Rev. D 37, 3406 (1988), URL http://link.aps.org/doi/10.1103/PhysRevD.37.3406.
- Copeland et al. (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D15, 1753 (2006), eprint hep-th/0603057.
- Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
- Carroll et al. (2004) S. M. Carroll, V. Duvvuri, M. Trodden, and M. S. Turner, Phys. Rev. D70, 043528 (2004), eprint astro-ph/0306438.
- Capozziello et al. (2003) S. Capozziello, S. Carloni, and A. Troisi, Recent Res. Dev. Astron. Astrophys. 1, 625 (2003), eprint astro-ph/0303041.
- Nojiri and Odintsov (2006) S. Nojiri and S. D. Odintsov, eConf C0602061, 06 (2006), [Int. J. Geom. Meth. Mod. Phys.4,115(2007)], eprint hep-th/0601213.
- Will (2006) C. M. Will, Living Rev. Rel. 9, 3 (2006), eprint gr-qc/0510072.
- Khoury (2010) J. Khoury (2010), eprint 1011.5909.
- Khoury (2013) J. Khoury (2013), eprint 1312.2006.
- Khoury and Weltman (2004a) J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004a), URL http://link.aps.org/doi/10.1103/PhysRevLett.93.171104.
- Khoury and Weltman (2004b) J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 (2004b), URL http://link.aps.org/doi/10.1103/PhysRevD.69.044026.
- Babichev et al. (2009) E. Babichev, C. Deffayet, and R. Ziour, Int. J. Mod. Phys. D18, 2147 (2009), eprint 0905.2943.
- Dvali et al. (2011) G. Dvali, G. F. Giudice, C. Gomez, and A. Kehagias, JHEP 08, 108 (2011), eprint 1010.1415.
- Vainshtein (1972) A. Vainshtein, Physics Letters B 39, 393 (1972), ISSN 0370-2693, URL http://www.sciencedirect.com/science/article/pii/0370269372901475.
- Hinterbichler and Khoury (2010) K. Hinterbichler and J. Khoury, Phys. Rev. Lett. 104, 231301 (2010), URL http://link.aps.org/doi/10.1103/PhysRevLett.104.231301.
- Olive and Pospelov (2008) K. A. Olive and M. Pospelov, Phys. Rev. D77, 043524 (2008), eprint 0709.3825.
- Abbott et al. (2005) T. Abbott et al. (Dark Energy Survey) (2005), eprint astro-ph/0510346.
- Comparat et al. (2012) J. Comparat, J.-P. Kneib, S. Escoffier, J. Zoubian, A. Ealet, et al., Mon.Not.Roy.Astron.Soc. 2012 (2012), eprint 1207.4321.
- Levi et al. (2013) M. Levi et al. (DESI collaboration) (2013), eprint 1308.0847.
- Abell et al. (2009) P. A. Abell et al. (LSST Science Collaborations, LSST Project) (2009), eprint 0912.0201.
- Laureijs et al. (2011) R. Laureijs et al. (EUCLID Collaboration) (2011), eprint 1110.3193.
- Spergel et al. (2013) D. Spergel, N. Gehrels, J. Breckinridge, M. Donahue, A. Dressler, et al. (2013), eprint 1305.5422.
- Zeldovich (1970) Ya. B. Zeldovich, Astron. Astrophys. 5, 84 (1970).
- Peebles (1980) P. J. E. Peebles, The large-scale structure of the universe (1980).
- Bouchet et al. (1995) F. R. Bouchet, S. Colombi, E. Hivon, and R. Juszkiewicz, Astron. Astrophys. 296, 575 (1995), eprint astro-ph/9406013.
- Tassev et al. (2013) S. Tassev, M. Zaldarriaga, and D. Eisenstein, JCAP 1306, 036 (2013), eprint 1301.0322.
- Brax et al. (2012) P. Brax, A.-C. Davis, B. Li, H. A. Winther, and G.-B. Zhao, JCAP 10, 002 (2012), eprint 1206.3568.
- Brax et al. (2013) P. Brax, A.-C. Davis, B. Li, H. A. Winther, and G.-B. Zhao, JCAP 1304, 029 (2013), eprint 1303.0007.
- Winther and Ferreira (2015) H. A. Winther and P. G. Ferreira, Phys. Rev. D91, 123507 (2015), eprint 1403.6492.
- Li and Barrow (2011) B. Li and J. D. Barrow, Phys. Rev. D83, 024007 (2011), eprint 1005.4231.
- Valkenburg and Hu (2015) W. Valkenburg and B. Hu, JCAP 1509, 054 (2015), eprint 1505.05865.
- Horndeski (1974) G. W. Horndeski, International Journal of Theoretical Physics 10, 363 (1974), ISSN 1572-9575, URL http://dx.doi.org/10.1007/BF01807638.
- Deffayet et al. (2011) C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, Phys. Rev. D 84, 064039 (2011), URL http://link.aps.org/doi/10.1103/PhysRevD.84.064039.
- Brax et al. (2012) P. Brax, A.-C. Davis, B. Li, and H. A. Winther, Phys. Rev. D86, 044015 (2012), eprint 1203.4812.
- Li and Zhao (2009) B. Li and H. Zhao, Phys. Rev. D 80, 044027 (2009), URL http://link.aps.org/doi/10.1103/PhysRevD.80.044027.
- Brax et al. (2008) P. Brax, C. van de Bruck, A.-C. Davis, and D. J. Shaw, Phys. Rev. D78, 104021 (2008), eprint 0806.3415.
- Hu and Sawicki (2007) W. Hu and I. Sawicki, Phys. Rev. D76, 064004 (2007), eprint 0705.1158.
- Zhao et al. (2011) G.-B. Zhao, B. Li, and K. Koyama, Phys. Rev. D83, 044007 (2011), eprint 1011.1257.
- Davis et al. (2012) A.-C. Davis, B. Li, D. F. Mota, and H. A. Winther, Astrophys. J. 748, 61 (2012), eprint 1108.3081.
- Klypin and Holtzman (1997) A. Klypin and J. Holtzman (1997), eprint astro-ph/9712217.
- Stabenau and Jain (2006) H. F. Stabenau and B. Jain, Phys. Rev. D74, 084007 (2006), eprint astro-ph/0604038.
- Laszlo and Bean (2008) I. Laszlo and R. Bean, Phys. Rev. D77, 024048 (2008), eprint 0709.0307.
- Khoury and Wyman (2009) J. Khoury and M. Wyman, Phys. Rev. D80, 064023 (2009), eprint 0903.1292.
- Lewis et al. (2000) A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), eprint astro-ph/9911177.
- Quinn et al. (1997) T. R. Quinn, N. Katz, J. Stadel, and G. Lake, Submitted to: Astrophys. J. (1997), eprint astro-ph/9710043.
- Scoccimarro (1998) R. Scoccimarro, Mon. Not. Roy. Astron. Soc. 299, 1097 (1998), eprint astro-ph/9711187.
- Smith et al. (2003) 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. Couchmann (VIRGO Consortium), Mon. Not. Roy. Astron. Soc. 341, 1311 (2003), eprint astro-ph/0207664.
- Kaiser (1987) N. Kaiser, Monthly Notices of the Royal Astronomical Society 227, 1 (1987), eprint http://mnras.oxfordjournals.org/content/227/1/1.full.pdf+html, URL http://mnras.oxfordjournals.org/content/227/1/1.abstract.
- Peacock and West (1992) J. A. Peacock and M. J. West, MNRAS 259, 494 (1992).
- Wang et al. (2013) Y. Wang, C.-H. Chuang, and C. M. Hirata, MNRAS 430, 2446 (2013), eprint 1211.0532.
- Davis and Peebles (1983) M. Davis and P. J. E. Peebles, Astrophys. J. 267, 465 (1983).
- Peacock and Dodds (1994) J. A. Peacock and S. J. Dodds, Mon. Not. Roy. Astron. Soc. 267, 1020 (1994), eprint astro-ph/9311057.
- Scoccimarro (2004) R. Scoccimarro, Phys. Rev. D70, 083007 (2004), eprint astro-ph/0407214.
- Angulo et al. (2008) R. Angulo, C. M. Baugh, C. S. Frenk, and C. G. Lacey, Mon. Not. Roy. Astron. Soc. 383, 755 (2008), eprint astro-ph/0702543.
- Kwan et al. (2012) J. Kwan, G. F. Lewis, and E. V. Linder, Astrophys. J. 748, 78 (2012), eprint 1105.1194.
- Jennings et al. (2012) E. Jennings, C. M. Baugh, B. Li, G.-B. Zhao, and K. Koyama, MNRAS 425, 2128 (2012), eprint 1205.2698.
- Behroozi et al. (2013) P. S. Behroozi, R. H. Wechsler, and H.-Y. Wu, Astrophys. J. 762, 109 (2013), eprint 1110.4372.
- Murray et al. (2013) S. Murray, C. Power, and A. Robotham (2013), eprint 1306.6721.
- Landy and Szalay (2002) S. D. Landy and A. S. Szalay, Astrophys. J. 579, 76 (2002), eprint astro-ph/0204499.
- Winther et al. (2017) H. A. Winther, K. Koyama, M. Manera, B. S. Wright, and G.-B. Zhao (2017), eprint arXiv:1703.00879.