Nonlinear structure formation in Nonlocal Gravity

# Nonlinear structure formation in Nonlocal Gravity

Alexandre Barreira Institute for Computational Cosmology, Department of Physics, Durham University, Durham DH1 3LE, U.K. Institute for Particle Physics Phenomenology, Department of Physics, Durham University, Durham DH1 3LE, U.K.    Baojiu Li Institute for Computational Cosmology, Department of Physics, Durham University, Durham DH1 3LE, U.K.    Wojciech A. Hellwing Institute for Computational Cosmology, Department of Physics, Durham University, Durham DH1 3LE, U.K. Interdisciplinary Centre for Mathematical and Computational Modeling (ICM), University of Warsaw, ul. Pawińskiego 5a, Warsaw, Poland    Carlton M. Baugh Institute for Computational Cosmology, Department of Physics, Durham University, Durham DH1 3LE, U.K.    Silvia Pascoli Institute for Particle Physics Phenomenology, Department of Physics, Durham University, Durham DH1 3LE, U.K.
###### Abstract

We study the nonlinear growth of structure in nonlocal gravity models with the aid of N-body simulation and the spherical collapse and halo models. We focus on a model in which the inverse-squared of the d’Alembertian operator acts on the Ricci scalar in the action. For fixed cosmological parameters, this model differs from by having a lower late-time expansion rate and an enhanced and time-dependent gravitational strength ( larger today). Compared to today, in the nonlocal model, massive haloes are slightly more abundant (by at ) and concentrated ( enhancement over a range of mass scales), but their linear bias remains almost unchanged. We find that the Sheth-Tormen formalism describes the mass function and halo bias very well, with little need for recalibration of free parameters. The fitting of the halo concentrations is however essential to ensure the good performance of the halo model on small scales. For , the amplitude of the nonlinear matter and velocity divergence power spectra exhibits a modest enhancement of to , compared to today. This suggests that this model might only be distinguishable from by future observational missions. We point out that the absence of a screening mechanism may lead to tensions with Solar System tests due to local time variations of the gravitational strength, although this is subject to assumptions about the local time evolution of background averaged quantities.

preprint: IPPP/14/ 70 DCPT/14/ 140

## I Introduction

It has been almost a century since Einstein proposed his theory of General Relativity (GR) which is still considered one of the main pillars of modern physics. The outstanding success of GR comes mostly from its ability to pass a number of stringent tests of gravity performed in the Solar System Will:2014kxa . When applied on cosmological scales, however, GR seems to lose some of its appeal as it requires the presence of some unknown form of dark energy in order to explain the observed accelerated expansion of the Universe. The simplest candidate for dark energy is a cosmological constant, , but the value of required to explain the observations lacks theoretical support. This has provided motivation for the proposal of alternative gravity models which attempt to reproduce cosmic acceleration without postulating the existence of dark energy. Furthermore, the fact that the laws of gravity have never been tested directly on scales larger than the Solar System justifies the exploration of such modifications to GR on cosmological scales. By understanding better the various types of observational signatures that different modified gravity models can leave on cosmological observables, one can improve the chance of identifying any departures from GR, or alternatively, extend the model’s observational success into a whole new regime. Currently, the study of modified gravity models is one of the most active areas of research in both theoretical and observational cosmology Clifton:2011jh ; Huterer:2013xky ; Jain:2013wgs ; Joyce:2014kja .

Here, we focus on a class of model that has attracted much attention recently, which is known as nonlocal gravity Woodard:2014iga . In these models, the modifications to gravity arise via the addition of nonlocal terms (i.e. which depend on more than one point in spacetime) to the Einstein field equations. These terms typically involve the inverse of the d’Alembertian operator, , acting on curvature tensors. To ensure causality, such terms must be defined with the aid of retarded Green functions (or propagators). However, it is well known that such retarded operators cannot be derived from standard action variational principles (see e.g. Sec. 2 of Ref. Woodard:2014iga for a discussion). One way around this is to specify the model in terms of its equations of motion and not in terms of its action. One may still consider a nonlocal action to derive a set of causal equations of motion, so long as in the end one replaces, by hand, all of the resulting operators by their retarded versions. Both of these approaches, however, imply that nonlocal models of gravity must be taken as purely phenomenological and should not be interpreted as fundamental theories. In general, one assumes that there is an unknown fundamental (local) quantum field theory of gravity, and the nonlocal model represents only an effective way of capturing the physics of the fundamental theory in some appropriate limit.

It was in the above spirit that Ref. Deser:2007jk proposed a popular nonlocal model of gravity capable of explaining cosmic acceleration. In this model, which has been extensively studied (see e.g. Refs. Deffayet:2009ca ; Deser:2013uya ; Woodard:2014iga ; Nojiri:2007uq ; Jhingan:2008ym ; Koivisto:2008xfa ; Koivisto:2008dh ; 2012PhRvD..85d4002E ; 2013CQGra..30c5002E ; 2013PhRvD..87b4003P ; Dodelson:2013sma and references therein), one adds the term to the Einstein-Hilbert action, where is the Ricci scalar and is a free function. As described in Ref. Woodard:2014iga , the function can be constructed in such a way that it takes on different values on the cosmological background and inside gravitationally bound systems. In particular, at the background level, can be tuned to reproduce -like expansion histories, but inside regions like the Solar System, one can assume that vanishes, thus recovering GR completely. This model, however, seems to run into tension with data sensitive to the growth rate of structure on large scales 2013PhRvD..87b4003P ; Dodelson:2013sma .

More recently, nonlocal terms have also been used to construct theories of massive gravity. An example of this is obtained by adding directly to the Einstein field equations a term like Maggiore:2013mea ; Foffa:2013vma ; Kehagias:2014sda ; Nesseris:2014mea , where is a mass scale and means the extraction of the transverse part (see also Refs. Jaccard:2013gla ; Foffa:2013sma ; Modesto:2013jea ; Ferreira:2013tqn for models in which acts on the Einstein and Ricci tensors). This model has no limit for the background evolution, but it can still match the current background expansion and growth rate of structure data with a similar goodness-of-fit Nesseris:2014mea . Furthermore, Ref. Kehagias:2014sda has investigated spherically symmetric static solutions in this model, concluding that it does not suffer from instabilities that usually plague theories of massive gravity. A similar model was proposed by Ref. Maggiore:2014sia , which is characterized by a term in the action (see Eq. (1)). Reference Dirian:2014ara showed that this model can reproduce current type Ia Supernovae (SNIa) data, although it also has no limit for the background expansion. The time evolution of linear matter density fluctuations in this model also differs from that in , but the work of Ref. Dirian:2014ara suggests that the differences between these two models are small enough to be only potentially distinguishable by future observational missions.

Here, we extend the previous work done for the model of Refs. Maggiore:2014sia ; Dirian:2014ara by examining its predictions in the nonlinear regime of structure formation. We achieve this by running a set of N-body simulations, which we use to analyse the model predictions for the nonlinear matter and velocity divergence power spectra, and also halo properties such as their abundance, bias and concentration. To the best of our knowledge this is the first time N-body simulations have been used to study the nonlinear regime of structure formation in nonlocal gravity cosmologies. N-body simulations are however computationally expensive to run. To overcome this, it is common to try to devise semi-analytical formulae that, motivated by simple physical assumptions, aim to reproduce the results from the simulations with the calibration of free parameters kept to a minimum. A popular example is given by the Sheth-Tormen halo mass function Sheth:1999mn ; Sheth:1999su ; Sheth:2001dp and its use in the halo model approach for the nonlinear matter power spectrum Cooray:2002dia . The performance of the halo model is well understood within , but less so in alternative gravity scenarios. In particular, the free parameters in these formulae might need substantial recalibration in models that differ significantly from (see e.g. Ref. 2014JCAP…04..029B ). One of our goals is to assess the performance of these analytical formulae in nonlocal gravity models.

This paper is organized as follows. In Sec. II we present the model and layout the equations relevant for the background evolution. We also derive the equations of motion for spherically symmetric configurations under the quasi-static and weak-field approximations. In Sec. III we present the formulae relevant for the calculation of the nonlinear matter power spectrum in the halo model formalism. In particular, we describe the Sheth-Tormen expressions for the halo mass function and bias, and define the Navarro-Frenk-White (NFW) halo concentration parameter. We also present the equations relevant for the linear evolution and spherical collapse of matter overdensities. Our results are presented in Sec. IV, where we discuss the results from the N-body simulations and compare them with the predictions from the analytical formulae. We also comment on the role that Solar System tests of gravity could play in setting the observational viability of this model. We summarize our findings in Sec. V.

In this paper, we work with the metric signature and use units in which the speed of light . Latin indices run over and Greek indices run over . We use interchangeably, where is the reduced Planck mass and is Newton’s constant.

## Ii The R□−2R nonlocal gravity model

### ii.1 Action and field equations

We consider the nonlocal gravity model of Refs. Maggiore:2014sia ; Dirian:2014ara , whose action is given by

 A=12κ∫dx4√−g[R−m26R□−2R−Lm], (1)

where is the determinant of the metric , is the Lagrangian density of the matter fluid, is the Ricci scalar and is the d’Alembertian operator. To facilitate the derivation of the field equations, and to solve them afterwards, it is convenient to introduce two auxiliary scalar fields defined as

 U = −□−1R, (2) S = −□−1U=□−2R. (3)

The solutions to Eqs. (2) and (3) can be obtained by evaluating the integrals

 U ≡ −□−1R = Uhom(x)−∫d4y√−g(y)G(x,y)R(y), S ≡ −□−1U = Shom(x)−∫d4y√−g(y)G(x,y)U(y),

where and are any solutions of the homogeneous equations and , respectively, and is any Green function of . The choice of the homogeneous solutions and of the Green function specify the meaning of the operator . To ensure causality, one should use the retarded version of the Green function, i.e., the solutions of (or ) should only be affected by the values of (or ) that lie in its past light-cone. The homogeneous solutions can be set to any value, which is typically zero, without any loss of generality. In principle, the model predictions can be obtained by solving Eqs. (II.1) and (II.1). However, it is convenient to use the fields and to cast the nonlocal action of Eq. (1) in the form of a local scalar-tensor theory Nojiri:2007uq ; Capozziello:2008gu ; Koshelev:2008ie as

 A =12κ∫dx4√−g[R−m26RS−ξ1(□U+R) (6) −ξ2(□S+U)−Lm],

where and are Lagrange multipliers. The field equations can then be written as

 Gμν−m26Kμν=κTμν, (7) □U=−R, (8) □S=−U, (9)

with

 Kμν≡2SGμν−2∇μ∇νS−2∇(μS∇ν)U +(2□S+∇αS∇αU−U22)gμν, (10)

and where is the energy-momentum tensor of the matter fluid. The use of the scalar fields and therefore allows one to obtain the solutions by solving a set of coupled differential equations, instead of the more intricate integral equations associated with the inversion of a differential operator. These two formulations are, however, not equivalent as explained with detail in many recent papers (see e.g. Refs.Koshelev:2008ie ; Koivisto:2009jn ; Barvinsky:2011rk ; Deser:2013uya ; Maggiore:2013mea ; Foffa:2013sma ; Foffa:2013vma ): Eqs. (7), (8) and (9) admit solutions that are not solutions of the original nonlocal problem. For instance, if is a solution of Eq. (8), then is also a solution for any , since (the same applies for the field and Eq. (9)). If one wishes the differential equations (7), (8) and (9) to describe the nonlocal model, then one must solve them with the one and only choice of initial conditions that is compatible with the choice of homogeneous solutions in Eqs. (II.1) and (II.1). All other initial conditions lead to spurious solutions and should not be considered.

### ii.2 Background equations

Throughout, we always work with the perturbed Friedmann-Roberston-Walker (FRW) line element in the Newtonian gauge,

 ds2=(1+2Ψ)dt2−a(t)2(1−2Φ)γijdxidxj, (11)

where is the cosmic scale factor ( is the redshift) and the gravitational potentials , are assumed to be functions of time and space. is the spatial sector of the metric, which is taken here to be flat.

At the level of the cosmological background (), the two Friedmann equations can be written as

 3H2 = κ¯ρm+κ¯ρde (12) −2˙H−3H2 = κ¯pm+κ¯pde, (13)

where we have encapsulated the effects of the nonlocal term into an effective background "dark energy" density, , and pressure , which are given, respectively, by

 κ¯ρde = m26[6¯SH2+6H˙¯S−˙¯U˙¯S−¯U22], (14) κ¯pde = −m26[2¯S(2˙H+3H2)+¨¯S +4H˙¯S+˙¯U˙¯S−¯U22].

Additionally, Eqs. (8) and (9) yield

 ¨¯U+3H˙¯U = 6(˙H+2H2), (16) ¨¯S+3H˙¯S = −¯U. (17)

In the above equations, a dot denotes a partial derivative w.r.t. physical time, , an overbar indicates that we are considering only the background average and is the Hubble expansion rate.

The background evolution in the model has to be obtained numerically. The differential equations are evolved starting from deep into the radiation dominated era () with initial conditions for the auxiliary fields . Note that, in the radiation era, the Ricci scalar vanishes (). Hence, from Eqs. (II.1) and (II.1) one sees that these initial conditions are indeed compatible with the choice . The value of the parameter is determined by a trial-and-error scheme to yield the value of that makes the Universe spatially flat, i.e., , where the subscripts , refer to radiation and matter, respectively, the subscript denotes present-day values, and is the present-day Hubble rate.

### ii.3 Spherically symmetric nonlinear equations

By assuming that the potentials and are spherically symmetric, one can write the and components of Eq. (7), and Eqs. (8) and (9), respectively, as

 2r2(r2Φ,r),r−m26[6SH2+4Sr2(r2Φ,r),r−2r2(r2S,r),r+2S,rΦ,r−S,rU,r−U22]=κ¯ρmδa2, (19) 2r(Φ,r−Ψ,r)−m26[4S˙Ha2+6SH2a2+4Sr(Φ,r−Ψ,r)+4S,rΦ,r−2S,rΨ,r−4S,rr+2S,rU,r−U22]=0, (20) 1r2(r2U,r),r+U,r(Ψ,r−Φ,r)=21r2(r2Ψ,r),r−41r2(r2Φ,r),r, (21) 1r2(r2S,r),r+S,r(Ψ,r−Φ,r)=U, (22)

where denotes a partial derivative w.r.t. the comoving radial coordinate . When writing Eqs. (19)-(22), we have already employed the following simplifying assumptions:

1. We have assumed the so-called quasi-static limit, under which one neglects the time derivatives of perturbed quantities, e.g., , where is the perturbed part of the auxiliary field;

2. We have also employed the so-called weak-field limit, which accounts for neglecting terms that involve and , and their first spatial derivatives, over those that involve their second spatial derivatives. For example, and .

The above equations still contain terms with and , because these terms contain the fields and , and up to now, we have not discussed the validity of applying these approximations to the auxiliary fields. However:

1. Equation (21) tells us that the field is of the same order as the scalar potentials, . Consequently, the above approximations also hold for ;

2. Equation (22) tells us that , which means we can also neglect all terms containing , and .

Under these considerations, the above equations simplify drastically. In particular, the only equation that remains relevant for the study of the spherical collapse of matter overdensities is Eq. (19), which can be written as:

 1r2(r2Φ,r),r=4πGeff¯ρmδa2, (23)

where

 Geff=G[1−m2¯S3]−1. (24)

Equation (23) is the same as in standard gravity, but with Newton’s constant replaced by the time-dependent gravitational strength, . This time dependence follows directly from the term in the field equations, Eq. (7), which in turn follows from the variation of the term in the action Eq. (6). The fact that depends only on time tells us that gravity is modified with equal strength everywhere, regardless of whether or not one is close to massive bodies or in high-density regions. This may bring into question the ability of this model to pass the stringent Solar System tests of gravity Will:2014kxa ; Babichev:2011iz ; 2012PhRvD..85b4023K . We come back to this discussion in Sec. IV.3. We note also that from Eq. (20), it follows that in the quasi-static and weak-field limits.

### ii.4 Model parameters

The results presented in this paper are for the cosmological parameter values listed in Table 1. These are the best-fitting parameters to a dataset that comprises the CMB data from the Planck satellite (both temperature and lensing) Ade:2013kta ; Ade:2013zuv ; Ade:2013tyw , and the BAO data from the 6df Beutler:2011hx , SDDS DR7 Padmanabhan:2012hf and BOSS DR9 Anderson:2012sa galaxy redshift surveys. The parameters were found by following the steps outlined in Ref. Barreira:2014jha , although in the latter, neutrino masses are also varied in the constraints. In this paper, however, we treat neutrinos as massless for simplicity.

The CMB temperature power spectra of the and models for the parameters listed in Table 1 are shown in Fig. 1. The model predictions were obtained with a suitably modified version of the CAMB code camb_notes . The derivation of the perturbed equations that enter the calculations in CAMB follows the steps of Ref. Barreira:2012kk , to which we refer the interested reader for details. The results in Fig. 1 shows that the model is able to fit the CMB data with a goodness-of-fit that is similar to that of . In fact, the model is in slightly better agreement with the data at low-, which is mostly determined by the Integrated Sachs-Wolfe (ISW) effect. However, the larger errorbars on these scales due to cosmic variance do not allow stringent constraints to be derived.

In Sec. IV we shall compare the results of the model with those of standard . In this paper, we are mostly interested in the phenomenology driven by the modifications to gravity in the model. This is why we shall use the same cosmological parameters for both models. A formal exploration of the constraints on the parameter space in the model is beyond the scope of the present paper (see Ref. swissgroup ).

## Iii Halo Model of the nonlinear matter power spectrum

In this section, we describe the halo model of the nonlinear matter power spectrum, as well as all of its ingredients. In particular, we define the halo mass function, linear halo bias and halo density profiles. We also present the equations that govern the linear growth and spherical collapse of structures.

### iii.1 Halo model

In the halo model approach, one assumes that all matter in the Universe lies within gravitationally bound structures (see Ref. Cooray:2002dia for a review). As a result, the two-point correlation function of the matter field can be decomposed into the contributions from the correlations between elements that lie in the same halo (the 1-halo term) and in different haloes (the 2-halo term). The power spectrum can also be decomposed in a similar way, and one can write

 Pk=P1hk+P2hk, (25)

where

 P1hk = ∫dMM¯ρ2m0dn(M)dlnM|u(k,M)|2, P2hk = I(k)2Pk,lin, (26)

are, respectively, the 1- and 2-halo terms, with

 I(k)=∫dM1¯ρm0dn(M)dlnMblin(M)|u(k,M)|. (27)

In Eqs.(25)-(27), is the present-day background (total) matter density; is the matter power spectrum obtained using linear theory; is the comoving wavenumber; is the mass function, which describes the comoving number density of haloes per differential logarithmic interval of mass; is the Fourier transform of the density profile of the haloes truncated at their size and normalized such that ; is the linear halo bias parameter. We model all these quantities in the remainder of this section, in which we follow the notation of Refs. Barreira:2013xea ; 2014JCAP…04..029B .

### iii.2 Halo mass function

We define the halo mass function as

 dn(M)dlnMdlnM=¯ρm0Mf(S)dS, (28)

where is the variance of the linear density field filtered on a comoving length scale ,

 S(R)≡σ2(R)=12π2∫k2Pk,lin~W2(k,R)dk. (29)

Here, is the Fourier transform of the filter, which we take as a top-hat in real space. The total mass enclosed by the filter is given by

 M=4π¯ρm0R3/3. (30)

In Eq. (28), describes the fraction of the total mass that resides in haloes whose variances lie within (or equivalently, whose masses lie within ) 111Note that the quantities , and can be related to one another via Eqs. (29) and (30). In this paper, we use these three quantities interchangeably when referring to the scale of the haloes.. Here, we use the Sheth-Tormen expression Sheth:1999mn ; Sheth:1999su ; Sheth:2001dp ,

 f(S)=A√q2πδcS3/2⎡⎣1+(qδ2cS)−p⎤⎦exp[−qδ2c2S],

where is a normalization constant fixed by the condition . is the critical initial overdensity for a spherical top-hat to collapse at redshift , extrapolated to with the CDM linear growth factor. This extrapolation is done purely to ensure that the resulting values of can be more readily compared to values in . Note that for consistency, in Eq. (29) is also the initial power spectrum of the specific model, evolved to with the CDM linear growth factor.

The Press-Schechter mass function 1974ApJ…187..425P is obtained by taking in Eq. (III.2). This choice is motivated by the spherical collapse model. However, Refs. Sheth:1999mn ; Sheth:1999su ; Sheth:2001dp showed that the choice of parameters (motivated by the ellipsoidal, instead of spherical collapse) leads to a better fit to the mass function measured from N-body simulations of CDM models (see also Ref. 2010ApJ…717..515M ). For alternative models, such as those with modified gravity, it is not necessarily true that the standard ST parameters () also provide a good description of the simulation results. For example, in Ref. 2014JCAP…04..029B , we demonstrated that the ST mass function can only provide a good fit to the simulation results of Galileon gravity models PhysRevD.79.064036 ; DeFelice:2010pv ; Deffayet:2009mn ; Barreira:2012kk after a recalibration of the parameters. In Sec. IV.4, we shall investigate the need for a similar recalibration in the model.

### iii.3 Linear halo bias

The linear halo bias parameter Mo:1995cs relates the clustering amplitude of haloes of mass to that of the total matter field on large length scales (),

 δhalo(M)=b(M)δmatter, (32)

where and are the density contrast of the distribution of haloes of mass and of the total matter field, respectively. On smaller length scales, where the matter overdensities become larger , Eq. (32) requires higher order corrections (see e.g. Fry:1992vr ).

Following the same derivation steps as in Ref. Barreira:2013xea , it is straightforward to show that the ST linear halo bias parameter can be written as

 b(M)=1+g(z)(qδ2c/S−1δc+2p/δc1+(qδ2c/S)p), (33)

with , where is the linear growth factor of a specific model defined as .

### iii.4 Halo density profiles

We adopt the NFW formula Navarro:1996gj to describe the radial density profile of dark matter haloes

 ρNFW(r)=ρsr/rs[1+r/rs]2, (34)

where and are often called the characteristic density and the scale radius of the halo. The mass of the NFW halo, , is obtained by integrating Eq. (34) up to some radius (the meaning of the subscript will become clear later)

 MΔ = ∫RΔ0dr4πr2ρNFW(r) = 4πρsR3Δc3Δ[ln(1+cΔ)−cΔ1+cΔ],

where is the concentration parameter.

In our simulations, we define the halo mass as

 MΔ=4π3Δ¯ρc0R3Δ, (36)

i.e., is the mass that lies inside a comoving radius , within which the mean density is times the critical density of the Universe at the present day, . In this paper, we take . By combining the two mass definitions of Eqs. (III.4) and (36), one gets as a function of :

 ρs=13Δ¯ρc0c3Δ[ln(1+cΔ)−cΔ1+cΔ]−1. (37)

The NFW profile then becomes fully specified by the values of , which are determined by direct fitting to the density profiles of the haloes from the simulations. Equivalently, and as is common practice in the literature, one can specify the concentration-mass relation , instead of . In the context of cosmologies, the (mean) concentration-mass relation is well fitted by a power law function over a certain mass range Bullock:1999he ; Neto:2007vq ; Maccio':2008xb ; Prada:2011jf . The same is true, for instance, for Galileon gravity models 2014JCAP…04..029B , although with very different fitting parameters. A proper assessment of the performance of the halo model prescription therefore requires us to fit the concentration-mass relation of the simulations as well. This is done in Sec. IV.6. Given the relation , then the NFW density profile becomes completely specified by the halo mass .

Finally, since it is the Fourier transform of the profiles, , and not the profiles themselves, that enter Eqs. (III.1) and (27), we simply mention that it is possible to show that

 uNFW(k,M)=∫RΔ0dr4πr2sinkrkrρNFW(r)MΔ (38) = 4πρsr3s{sin(krs)M[Si([1+cΔ]krs)−Si(krs)] +cos(krs)M[Ci([1+cΔ]krs)−Ci(krs)] −sin(cΔkrs)M(1+cΔ)krs},

where and . Note that , as required by its normalization.

### iii.5 Linear growth factor and spherical collapse dynamics

In order to use the ST formulae for the mass function and linear halo bias one still needs to specify and solve the equations that determine the threshold density and the evolution of the linear overdensities. For scales well within the horizon, the linear (small) density contrast is governed by

 ¨δlin+2H˙δlin−4πGeff(a)¯ρmδlin=0, (39)

or equivalently, by changing the time variable to , by

 D′′+(E′E+2)D′−32Geff(a)GΩm0e−3NE2 = 0, (40)

where we have used that and a prime denotes a derivative w.r.t. . The initial conditions are set up at using the known matter dominated solution . The model changes the way structure grows on large scales compared to via its modifications to and .

We have defined as the linearly extrapolated value (using the CDM linear growth factor) of the initial overdensity of a spherical region for it to collapse at a given redshift, . To determine , we consider the evolution equation of the physical radius of the spherical halo at time , which satisfies the Euler equation

 ¨ζζ−(˙H+H2) = −Φ,ζζ=−Geff(a)GH20Ωm0δa−32, (41)

where the last equality follows from integrating Eq. (23) over . Changing the time variable to and defining , Eq. (41) becomes

 y′′ + (E′E+2)y′ (42) + Geff(a)GΩm0e−3N2E2(y−3−1)y=0,

where we have used , which follows from mass conservation 222Explicitly, one has .. The initial conditions are set up as and (here, is the linear density contrast at the initial time). The value of is then determined by finding the value of the initial density that leads to collapse (, ) at redshift , evolving this afterwards until today using the CDM linear growth factor.

As we have noted above, in the model, the modifications to gravity are time dependent only. In other words, the value of is the same on large and on small scales. This is different, for instance, from the case of Galileon gravity. In the latter, the nonlinearities of the Vainshtein screening mechanism suppress the effective gravitational strength felt by a test particle that lies within a certain radius (known as Vainshtein radius) from a matter source. In Eq. (42), this would be simply taken into account by replacing with Barreira:2013xea . The picture becomes more complicated in the case of modified gravity models which employ chameleon-type screening mechanisms Khoury:2003aq ; Khoury:2003rn . In these models, also depends on the size (or mass) of the halo and on the gravitational potentials in the environment where the halo forms. This requires a generalization of the spherical collapse formalism for these models, which has been developed by Ref. Li:2011qda .

## Iv Results

### iv.1 N-body simulations summary

Our simulations were performed with a modified version of the publicly available RAMSES N-body code Teyssier:2001cp . RAMSES is an Adaptive Mesh Refinement (AMR) code, which solves the Poisson equation on a grid that refines itself when the effective number of particles within a given grid cell exceeds a user-specified threshold, . Our modifications to the code consist of (i) changing the routines that compute the background expansion rate to interpolate the model expansion rate from a pre-computed table generated elsewhere; (ii) re-scaling the total force felt by the particles in the simulation by , whose values are also interpolated from a table generated beforehand.

In the following sections we show the N-body simulation results obtained for three models. We simulate the "full" model of action Eq. (1), whose expansion history and are given by Eqs. (12) and (24), respectively. We also simulate a standard model and a model with the same expansion history as the model, but with . We call the latter model , and comparing its results to allows us to pinpoint the impact of the modified alone on the growth of structure. The specific impact of the modified can then be measured by comparing the results from the "full" model simulations with those from . Table 2 summarizes the models we consider in this paper.

We simulate all models on a cubic box of size with dark matter particles. We take as the grid refinement criterion. The initial conditions are set up at , using the linear matter power spectrum with the parameters of Table 1. For each model, we simulate five realizations of the initial conditions (generated using different random seeds), which we use to construct errorbars for the simulation results by determining the variance across the realizations.

Finally, we simply note that the modifications to RAMSES needed to simulate the model are trivial compared to those that are necessary to simulate models such as Li:2011vk ; Puchwein:2013lza ; Llinares:2013jza ; Oyaizu:2008sr ; 2011PhRvD..83d4007Z ; Hellwing2013b , Symmetron 2012ApJ…748…61D ; Llinares:2013jua or Galileons 2009PhRvD..80j4005C ; Schmidt:2009sg ; Li:2013nua ; Wyman:2013jaa ; Barreira:2013eea ; Li:2013tda . In these, because of the screening mechanisms (which introduce density and scale dependencies of the total force), additional solvers are needed for the (nonlinear) equations of the extra scalar degrees of freedom. In the ECOSMOG code Li:2011vk (also based on RAMSES), these equations are solved via Gauss-Seidel relaxations on the AMR grid, which makes the simulations significantly more time consuming. We note also that recently, Ref. Winther:2014cia has proposed a new and faster scheme to simulate screened modified gravity in the mildly non-linear regime. This scheme uses the linear theory result, but combines it with a screening factor computed analytically assuming spherical symmetry, which helps speed up the calculations without sacrificing the accuracy on mildly nonlinear scales too much.

### iv.2 Linear growth and δc curves

Before discussing the results from the simulations, it is instructive to look at the model predictions for the linear growth rate of structure and for the time dependence of the critical density .

From top to bottom, Fig. 2 shows the time evolution of the fractional difference of the expansion rate relative to , , the effective gravitational strength and the fractional difference of the squared linear density contrast relative to , . The expansion rate in the model is lower than in for . This reduces the amount of Hubble friction and therefore boosts the linear growth rate. The gravitational strength in the model starts growing after , being approximately larger than in GR at the present day. This also boosts the linear growth of structure, but has a smaller impact compared to the effect of the lower expansion rate. This is seen by noting that the differences between and in the bottom panel are larger than the differences between and the model.

Figure 3 shows the time dependence of . In the top panel, all models exhibit the standard result that decreases with time, i.e., the initial overdensity of the spherical top-hat should be smaller, if the collapse is to occur at later times. Compared to , at late times (), the and models predict lower values for . This is as expected since structure formation is boosted at late times in these models, and as a result, this needs to be compensated by smaller values of the initial overdensities for the collapse to occur at the same epoch as in . Just like in the case of the linear growth rate, the differences w.r.t. the results are mainly affected by the lower expansion rate, and not by the larger values of . At earlier times (), all models have essentially the same expansion rate and gravitational strength, and as a result, the values of are roughly the same. Table 3 shows the values of at , and , for the three models.

### iv.3 Interpretation of the constraints from Solar System tests of gravity

The absence of a screening mechanism in the model may raise concerns about the ability of the model to satisfy Solar System constraints Will:2014kxa . For instance, for the parameters of Table 1, the model predicts that the rate of change of the gravitational strength today, , is

 ˙GeffG=H0ddN(GeffG)≈92×10−13 yrs−1, (43)

which is at odds with the observational contraint , obtained from Lunar Laser Ranging experiments Williams:2004qba . Hence, it seems that this type of local constraints can play a crucial role in determining the observational viability of the model, potentially ruling it out (see e.g. Refs. Babichev:2011iz ; 2012PhRvD..85b4023K for a similar conclusion, but in the context of other models).

It is interesting to contrast this result with that of the nonlocal model of Ref. Deser:2007jk , which we call here the model (for brevity), where is a free function that appears in the action and . The equations of motion of this model can be schematically written as

 Gμν[1+χ(X)]+ΔGμν=Tμν, (44)

where encapsulates all the extra terms that are not proportional to and the factor is given by

 χ=f(X)+□−1[RdfdX(X)]. (45)

For the purpose of our discussion, it is sufficient to look only at the effect of in Eq. (44). This rescales the gravitational strength as

 GeffG={1+χ}−1, (46)

which is similar to the effect of in the model. There is, however, one very important difference associated with the fact that in the case of the model, one has the freedom to choose the functional form of the terms that rescale . To be explicit, we write the argument of as

 X=□−1R=□−1¯R+□−1δR, (47)

where and are, respectively, the background and spatially perturbed part of . As explained in Ref. Woodard:2014iga , the relative size of and is different in different regimes. At the background level, and so the operator acts only on . On the other hand, within gravitationally bound objects we have . Now recall that the covariant operator acts with different signs on purely time- and space-dependent quantities 333For instance, in flat four-dimensional Minkowski space we have .. As a result, the sign of on the background differs from that within bound systems, such as galaxies or our Solar System. This can be exploited to tune the function in such a way that it vanishes when the sign of is that which corresponds to bound systems. In this way, and one recovers GR completely 444In Eq. (44), also vanishes if .. When takes the sign that corresponds to the background, then the function is tuned to reproduce a desired expansion history, typically . In the case of the model, is fixed to be and one does not have the freedom to set it to zero inside bound objects. Consequently, the time-dependent part of is always present in Eq. (24), which could potentially lead to a time-dependent gravitational strength that is at odds with the current constraints.

For completeness, one should be aware of a caveat. In the above reasoning, we have always assumed that the line element of Eq. (11) is a good description of the geometry of the Solar System. The question here is whether or not the factor should be included in the spatial sector of the metric when describing the Solar System. This is crucial as the presence of in Eq. (11) determines if varies with time or not. If is considered, then varies with time and is time-varying as well. In this way, the model fails the Solar System tests. On the other hand, if one does not consider in the metric, then is forcibly constant, and there are no apparent observational tensions. Such a static analysis was indeed performed by Refs. Kehagias:2014sda ; Maggiore:2014sia , where it was shown that the model can cope well with the local constraints.

This boils down to determining the impact of the global expansion of the Universe on local scales. It is not clear to us that if a field is varying on a time-evolving background, then it should not do so in a small perturbation around that background. However, we acknowledge this is an open question to address, and such study is beyond the scope of the present paper. In what follows, we limit ourselves to assuming that Eq. (24) holds on all scales, but focus only on the cosmological (rather than local) interpretation of the results.

### iv.4 Halo mass function

Our results for the cumulative mass function of the (black), (red) and (blue) models are shown in Fig. 4 at , and . The symbols show the simulation results obtained with the halo catalogues we built using the Rockstar halo finder Behroozi:2011ju . The results in the figure correspond to catalogues with subhaloes filtered out. The lines show the ST analytical prediction (Eqs. (28) and (III.2)) computed for the fitted (solid lines) and standard (dashed lines) ST parameters of Table 4. The parameters were fitted for all of the epochs shown, using the corresponding values of in Table 3. From the figure, one notes that although performing the fitting helps to improve the accuracy of the analytical formulae, overall the use of the standard values for provides a fair estimate of the halo abundances in the model, and of its relative difference w.r.t. . This is not the case, for instance, in Galileon gravity models, for which Ref. 2014JCAP…04..029B has found that it is necessary to recalibrate substantially the values of if the ST mass function is to provide a reasonable estimate of the effects of the modifications to gravity. In the case of the model, the fact that the standard values of work reasonably well means that the modifications in the model, relative to , are mild enough for its effects on the mass function to be well captured by the differences in .

At , the mass function of the model shows an enhancement at the high-mass end (), and a suppression at the low-mass end (), relative to . This is what one expects in hierarchical models of structure formation if the growth rate of structure is boosted, as smaller mass objects are assembled more efficiently to form larger structures, leaving fewer of them. The effects of the enhanced maintain this qualitative picture, but change it quantitatively. More explicitly, the mass scale below which the mass function drops below that of is smaller than the mass range probed by our simulations; and the enhancement of the number density of massive haloes is more pronounced. In particular, compared to , haloes with masses are