Constraining chameleon models with cosmology
aa r X i v : . [ a s t r o - ph . C O ] J a n Constraining chameleon models with cosmology
Lucas Lombriser Institute for Astronomy, University of Edinburgh,Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, U.K. (Dated: June 4, 2018)Chameleon fields may modify gravity on cluster scales while recovering general relativity locally.This article reviews signatures of chameleon modifications in the nonlinear cosmological structure,comparing different techniques to model them, summarising the current state of observational con-straints, and concluding with an outlook on prospective constraints from future observations andapplications of the analytic tools developed in the process to more general scalar-tensor theories.Particular focus is given to the Hu-Sawicki and designer models of f ( R ) gravity. I. INTRODUCTION
Attempts to unify general relativity with the standardmodel interactions typically introduce an effective scalarfield in the low-energy limit in addition to the gravita-tional tensor field. This fifth element may couple mini-mally or nonminimally to the matter fields and can actas an alternative to the cosmological constant as expla-nation for the observed late-time acceleration of our Uni-verse. A fifth force originating from the nonminimal cou-pling consequently modifies the gravitational interactionsbetween the matter fields. Local observations place tightconstraints on the existence of a fifth force [1]. However,these constraints are alleviated for scalar field potentialsthat yield a scalar field inversely dependent on curva-ture such as is realised in chameleon models [2–6]. Inthis case, the extra force is suppressed in high-densityregions like the Solar System but at the cost of returningcosmic acceleration to be driven by the contribution of adark energy component or a cosmological constant ratherthan originating from a fundamental modified gravity ef-fect of the model [7]. Meanwhile, gravitational forces re-main enhanced at low curvature and below the Comptonwavelength of the scalar field, causing an increase in thegrowth of structure and rendering nonlinear cosmologicalstructures a vital regime for testing gravity.This article reviews the signatures of chameleon fieldsin the cosmological small-scale structure, summarisingdifferent techniques to model the formation of thesestructures and current constraints on the chameleonfield amplitude and coupling strength from observations.Thereby, scalar-tensor models are considered that havea constant Brans-Dicke parameter ω , match the ΛCDMbackground expansion history, and exhibit a chameleonsuppression mechanism of the enhanced gravitationalforce in high-density regions. A special emphasis is givento the Hu-Sawicki [8] and designer [9] f ( R ) models [10–15], to which the scalar-tensor model can be reduced inthe case of ω = 0. These models have been particu-larly well studied with constraints from current and ex-pected from future observations reported, for instance,in Refs. [8, 9, 16–81]. The currently strongest boundson the chameleon modifications are inferred from com-paring nearby distance measurements in a sample of un- screened dwarf galaxies [50] as well as from the transitionof the scalar field required in the galactic halo to interpo-late between the high curvature regime of the Solar Sys-tem, where the chameleon mechanism suppresses forcemodifications, and the low curvature of the large-scalestructure, where gravity is modified [8, 75]. Cosmolog-ical observables such as cluster profiles [46] and abun-dance [29, 32, 38], galaxy power spectra [79], redshift-space distortions [33, 54], or the combination of gas andweak lensing measurements of a cluster [77] can be usedto place independent and strong bounds on the gravita-tional modifications.Local and astrophysical probes yield constraints thatare 2-3 orders of magnitude stronger than what is inferredfrom cosmological observations. It is worth emphasising,however, that they test the coupling of the scalar field tobaryons, whereas cosmological probes typically rely on acoupling of the scalar field to dark matter only, except forthe constraints inferred using the cluster gas in Ref. [77].The separation of coupling strengths between the differ-ent matter components can be made explicit in the Ein-stein frame, in which case cosmological constraints canbe regarded as independent of the local and astrophys-ical bounds. Furthermore, dark chameleon fields thatonly couple to dark matter may alleviate problems aris-ing through quantum particle production [82, 83] due tohigh-energy fluctuations in the early universe or largequantum corrections to the scalar field potential in lab-oratory environments [84].Hence, cosmological scales remain an interestingregime for constraining potential scalar field contribu-tions. The chameleon mechanism is a nonlinear effect,however, that complicates the description of the cos-mological small-scale structure. N -body simulations ofchameleon models provide an essential tool for study-ing the nonlinear regime of structure formation andthe suppression of gravitational modifications [85–92].These simulations are, however, computationally consid-erably more expensive than in the Newtonian scenario.Hence, for the comparison of theory to observations, al-lowing a full exploration of the cosmological parameterspace involved, the development of more efficient mod-elling techniques for the cosmological small-scale struc-ture is a necessity. Different approaches have been pro-posed based on phenomenological frameworks and fittingfunctions to simulations [93–96], analytical and numer-ical approximations [48, 77, 97–99], the spherical col-lapse model [75, 100–103], excursion set theory [102–106], the halo model [46, 75, 100], and perturbation the-ory [107, 108]. These different methods shall be sum-marised and compared here.Sec. II reviews chameleon gravity in the context ofmore general scalar-tensor theories and discusses the Hu-Sawicki and designer f ( R ) models. Sec. III is devoted tothe formation of large-scale structure in chameleon mod-els and its description using linear cosmological perturba-tion theory in the quasistatic limit, dark matter N -bodysimulations, and the spherical collapse model. It fur-thermore discusses different modelling techniques for thematter power spectrum and the properties of chameleonclusters such as the halo mass function and linear halobias as well as the cluster profiles of the matter density,scalar field, and dynamical mass. In Sec. IV, current con-straints on chameleon models, in specific on f ( R ) grav-ity and from cosmological observations, are summarised.Finally, Sec. V provides an outlook and discussion ofprospective constraints from future observations and ap-plications of the modelling techniques developed for thechameleon modifications to more general scalar-tensortheories, before Sec. VI concludes the review. II. CHAMELEON MODELS
The modified gravity models studied here are bestviewed as a special limit of the more general scalar-tensortheory defined by the Horndeski extension [109, 110] tothe Einstein-Hilbert action S = 12 κ Z d x √− g (cid:26) G ( ϕ, X ) − G ( ϕ, X ) (cid:3) ϕ + G ( ϕ, X ) R + ∂G ∂X (cid:2) ( (cid:3) ϕ ) − ( ∇ µ ∇ ν ϕ ) (cid:3) + G ( ϕ, X ) G µν ∇ µ ∇ ν ϕ − ∂G ∂X (cid:2) ( (cid:3) ϕ ) − (cid:3) ϕ ( ∇ µ ∇ ν ϕ ) + 2( ∇ µ ∇ ν ϕ ) (cid:3)(cid:27) + S m [ ψ m ; g µν ] , (1)where X ≡ − ( ∂ µ ϕ ) , R is the Ricci scalar, G µν is theEinstein tensor, S m is the matter action with matterfields ψ m , κ ≡ π G with the bare gravitational coupling G , and natural units are assumed here and throughoutthe article. The scalar field ϕ is coupled to the metric g µν via the covariant derivatives, R , and G µν . Eq. (1) de-fines the most general scalar-tensor theory for which theEuler-Lagrange equations involve at most second orderderivatives of the fields. Hereby, G − are free functionsof ϕ and X .In the following, we specialise to a subclass of the Horn-deski action, the Jordan-Brans-Dicke models, which are embedded in Eq. (1) via the definitions G ( ϕ, X ) ≡ − (cid:20) U ( ϕ ) − ω ( ϕ ) ϕ X (cid:21) ,G ( ϕ, X ) ≡ G ( ϕ, X ) ≡ ,G ( ϕ, X ) ≡ ϕ, (2)where U ( ϕ ) is the scalar field potential and ω ( ϕ ) is theBrans-Dicke parameter determining the kinetic coupling,assumed to be constant in the following with ω > − / G − terms, and the symmetron mod-els [112], which can be represented by appropriate choicesof ω ( ϕ ) and U ( ϕ ).The Jordan frame action defined by Eqs. (1) and (2)can be recast in the Einstein frame using the transforma-tions ˜ g µν ≡ ϕ g µν , (3) (cid:18) d φ d ϕ (cid:19) ≡ κ ωϕ , (4) A ( φ ) ≡ ϕ − / , (5) V ( φ ) ≡ U ( ϕ ) κ ϕ , (6)such that S = Z d x p − ˜ g " ˜ R κ − ∂ µ φ ∂ µ φ − V ( φ ) + S m (cid:2) ψ m ; A ( φ )˜ g µν (cid:3) , (7)where here and throughout the article, tildes denotequantities in the Einstein frame.It is worth emphasising that if considering the Einsteinframe action Eq. (7) to be the fundamental action defin-ing the chameleon model, the fifth element may not berestricted to couple to the different matter componentswith the same coupling strength, for instance, discrimi-nating between the baryonic components (b) and the colddark matter (c). In this case, the matter action in Eq. (7)is replaced by S m (cid:2) ψ m ; A ( φ )˜ g µν (cid:3) → S b+c with S b+c = S b (cid:2) ψ b ; A ( φ )˜ g µν (cid:3) + S c (cid:2) ψ c ; A ( φ )˜ g µν (cid:3) , (8)where one can define the coupling strengths β i to the dif-ferent matter components as A b ( φ ) = exp( β b κ φ ) and A c ( φ ) = exp( β c κ φ ). This underlines the importanceof cosmological tests of scalar field couplings as comple-ment to the local tests. For instance, in the scenario inwhich β b = 0, the scalar field is minimally coupled tothe baryons and Solar System or astrophysical tests asin Refs. [8, 50, 75] do not constrain β c = β , whereas themodel can still be constrained using cosmological obser-vations which test the distribution of the cold dark mat-ter. This review restricts to models in which the couplingstrengths to the different matter components are assumedequal, corresponding to β b = β c = β with constant β and β κ φ ≪ φ = 1 κ r ω ϕ + φ , (9)where in the following φ = 0. Variation of the actionEq. (7) with respect to φ yields the scalar field equation˜ (cid:3) φ = κ √ ω ˜ T + V ′ ( φ ) ≡ V ′ eff ( φ ) , (10)where V eff ( φ ) is an effective potential governing the dy-namics of φ and the energy-momentum tensor is given by˜ T = A ( φ ) T . For a scalar field with ϕ ≃ V ′ eff ( φ ) = 0, Eq. (10) becomesdd ϕ U ( ϕ ) ≃
12 ( κ ρ m + 4 U ) ≃ ˜ R ≃ R , (11)assuming dominance of the matter energy density ρ m andfurther requiring ( κ ρ m + 4 U ) ≫ (3 + 2 ω )( ∂ µ ϕ ) / A. Chameleon mechanism
With the appropriate choice of U ( ϕ ), the model speci-fied by Eqs. (1) and (2) produces cosmic acceleration andintroduces a modification of gravity, which can be sup-pressed in high-density regions in order to satisfy SolarSystem constraints. These demands on U ( ϕ ) can be sat-isfied with the ansatz U = Λ + U α (1 − ϕ ) α , where α isa positive constant such that in the limit of ϕ = 1, themodel reduces to ΛCDM. Hereby, Λ may be viewed as aneffective cosmological constant, valid as description in thelimit of large enough curvature of a more general U ( ϕ ),as is the case in the Hu-Sawicki f ( R ) models discussedin Sec. II B. On the other hand, chameleon models withequivalent coupling to the different matter fields cannotaccommodate cosmic acceleration as a genuine modifiedgravity effect while simultaneously satisfying Solar Sys-tem constraints [7]. Self-acceleration should be intro-duced as a consequence of the transformation Eq. (3)from a non-accelerating expansion in Einstein frame tothe Jordan frame, which is not the case here and cosmicacceleration is driven by a cosmological constant as in ouransatz for U ( ϕ ) or a dark energy contribution as, for in-stance, in the designer f ( R ) model discussed in Sec. II B.With κ ρ m ≫ −√ ω ˜ ∇ φ in high-density regions,assuming the quasistatic limit, the scalar field is in theminimum of the effective potential V eff ( φ ) in Eq. (10),and hence, U ϕ ≃ R/ α =1. Thus, in this case, and with the effective potential alsominimised in the background when | ω j ¯ ϕ ( i ) | ≪ j = 0 , i = 1 , ϕ ≃ ϕ − (cid:18) ¯ R R (cid:19) / (1 − α ) , (12) U ( ϕ ) ≃ Λ − ¯ R α (1 − ϕ ) α (1 − ¯ ϕ ) α − , (13)where subscripts of zero and overbars denote quantitiesevaluated at present time, a ≡
1, and in the background,respectively, here and throughout this review. Thus, inthe high-curvature regime, R ≫ ¯ R , for α <
1, the scalarfield is suppressed, ( ϕ − ≃ κ/ √ ω φ ≃ U ≃ Λ, and hence, reproduce a ΛCDMexpansion history with ∆ H ∼ O (1 − ¯ ϕ ), one further re-quires α ≫ | ¯ ϕ − | and | ω j || − ¯ ϕ | ≪ (1 − α ) such that | ω j ϕ ( i ) | ≪ | ¯ ϕ − | .
56 + 4 ω × − (14)with very weak dependence on α , which, however, is con-strained by the requirement that ¯ U ≈ Λ, implying that α ≫ − (6 + 4 ω ) − .Whereas the chameleon screening mechanism relies onthe scalar field mass becoming large in high-density re-gions, requiring the chameleon model to be valid as aneffective field theory with one-loop quantum correctionsto V ( φ ) not exceeding the classical scalar field poten-tial, places an upper bound on the mass of the field ofthe order of O (10 − ) eV within a laboratory environ-ment [84]. The corresponding minimal scalar field massset by the Solar System constraint in Eq. (14), however,is & O (1 −
10) eV. Hence, the effective field theory in-terpretation of the models considered in Eq. (13) is notquantum-stable on these scales. On the other hand, theclassical chameleon models of Eq. (13) have no effect inE¨ot-Wash type laboratory experiments [113].Besides the Solar System bounds, chameleon modelscan also be tightly constrained by astrophysical and cos-mological probes. Hereby, it is important to note thatin addition to the suppression of modifications in high-density regions due to the chameleon mechanism, onscales larger than the background Compton wavelengthtoday, ¯ m − ∼ r ω − α − ¯ ϕ − Mpc (15)(see Sec. III A), extra forces become Yukawa suppressed.Thus, with the local constraints in Eq. (14), this impliesthat ¯ m − . (1 − α ) − / Mpc and that modified gravityeffects are limited to nonlinear cosmological structures(cf. [7]).Among the strongest constraints on ¯ ϕ are the boundsthat can be inferred from astrophysical observations suchas from analysing different distance indicators in un-screened dwarf galaxies [50] or cosmological tests such asfrom the comparison of gas to weak lensing measurementsin the Coma cluster [77]. Fig. 1 summarises and com-pares these constraints from the different astronomicalregimes, including the local bounds in Eq. (14), as func-tion of the coupling strength ω and present backgroundfield amplitude ¯ ϕ . The corresponding constraints in thelimit of f ( R ) gravity (see Sec. II B) are indicated by thevertical dotted line. Note that these results depend onthe requirement that β b = β c and do not apply to scenar-ios where β b = 0, in which case, however, modificationscan, for instance, be constrained through signatures inthe observed cluster abundance [29, 32, 38, 46]. B. f ( R ) gravity In f ( R ) gravity, the Einstein-Hilbert action is sup-plemented with a free nonlinear function of the Ricciscalar, replacing R → R + f ( R ) in the Lagrangian den-sity [10–15]. This modification can be embedded in thethe Jordan-Brans-Dicke subclass of the Horndeski actiondefined by Eqs. (1) and (2) with ω ≡ ϕ ≡ f / d R ≡ f R , and scalar field potential U ≡ ( Rf R − f ) /
2. Twoparticularly well studied classes of f ( R ) models are theHu-Sawicki and designer forms.The Hu-Sawicki model [8] is described by the func-tional form f ( R ) = − H c (cid:0) R/H (cid:1) n c ( R/H ) n + 1 (16)with H ≡ κ ¯ ρ m0 /
3. The parameters c , c , and n are the degrees of freedom of the model, which can beset in order to match the ΛCDM background expansionhistory and satisfy Solar System constraints through thechameleon suppression mechanism. For sufficiently largecurvature, c /n R ≫ H , Eq. (16) simplifies to f ( R ) = − c c H − f R n ¯ R n +10 R n (17)with f R ≡ f R ( ¯ R ). Requiring equivalence with ΛCDMwhen | f R | → c c H = 2 κ ¯ ρ Λ . (18)Note that the scalar-tensor model defined by the combi-nation of Eqs. (1) and (2) with the scalar field potentialEq. (13) reduces to the Hu-Sawicki f ( R ) gravity model,Eq. (17), in the limit of ω = 0, where α = n/ ( n + 1).The corresponding normalised scalaron field f R /f R isshown as a function of R/ ¯ R for different values of α inFig. 1. The scalar field mass in the background evaluatedtoday ¯ m can be related to | f R | following Eq. (25) with¯ m ≃ − (4 − m )(1 − α ) | f R | − h Mpc − . Another well studied class of f ( R ) models are the de-signer models [9, 114–117]. In this case, f ( R ) is recon-structed from a predefined background expansion his-tory, which here, shall be given by the matter-dominatedΛCDM Hubble parameter H = κ (¯ ρ m + ¯ ρ Λ ) /
3. Us-ing this requirement in the Friedmann equations of f ( R )gravity yields an inhomogeneous second-order differentialequation for f ( R ), f ′′ − (cid:20) H ′′ H + ¯ R ′′ ¯ R ′ (cid:21) f ′ + ¯ R ′ H f = − H (1 − Ω m ) ¯ R ′ H , (19)where if not otherwise specified, primes denote deriva-tives with respect to ln a here and throughout the arti-cle. Eq. (19) can be solved numerically with the initialconditions f (ln a i ) = A H a p i − H (1 − Ω m ) , (20) f ′ (ln a i ) = p A H a p i , (21)where p = ( − √ /
4, set at an initial time a i ≪ A then characterises a particular solution inthe set of functions f ( R ) that exactly recover a ΛCDMbackground expansion history. Rather than characteris-ing the solutions by an initial condition A at an arbitraryredshift, they can be defined by f R as in the Hu-Sawickimodel or by the Compton wavelength parameter [9] B = f RR ( ¯ R )1 + f R ( ¯ R ) ¯ R ′ HH ′ (22)at its present value B ≡ B ( a = 1).Fig. 1 compares the normalised scalaron field f R /f R of the designer model, using the cosmological parametersdefined in Sec. III B, with the normalised scalaron fieldof the Hu-Sawicki model. From Fig. 1 it is clear that al-though both models reproduce the ΛCDM background,one exactly, the other one approximately, the scalar fieldsdo not match beyond ¯ R and the modifications of grav-ity are similar but not equivalent. The designer modelcannot be described by the Hu-Sawicki model using aconstant value of α . Note that in the limit of B = 0 orequivalently f R = 0, both the designer and Hu-Sawicki f ( R ) scenarios reduce to the ΛCDM model, i.e., includ-ing the perturbations to the background. III. CHAMELEON COSMOLOGY
Whereas the chameleon models discussed in Sec. IImatch the ΛCDM background expansion history and re-cover general relativity in the Solar System, the cos-mological structure on scales of a few megaparsecs andless remains modified. Importantly, this restriction ofgravitational modifications to small-scale cosmology isinferred from the local bounds on the coupling of the - - - - Ω+ (cid:144) - j f H R L CosmologySolar System Astrophysics - - - (cid:144) R f R (cid:144) f R Designer modelHu - Sawicki model
Α =
Α =
Α =
Α =
FIG. 1.
Left panel:
Local [75], astrophysical [50], and cosmological [77] constraints on the present amplitude of the backgroundchameleon field ¯ ϕ , depending on the Brans-Dicke parameter ω , which is assumed constant and describes the coupling of thechameleon field to matter. The vertical dotted line indicates the limit of f ( R ) gravity. Right panel:
Comparison of theHu-Sawicki and designer f ( R ) model scalaron fields. scalar field to baryons. Cosmological probes, however,typically test the distribution of dark matter. Hence,besides providing independent constraints on the samegravitational model at different scales of our Universe,they may also serve as incompensable tests of the pres-ence of scalar fields that couple nonminimally to the darkmatter but may not couple equally strong to the bary-onic components. Thus, it is important to study thesignatures of chameleon fields on the large-scale struc-ture even if the amplitudes and coupling strengths forthe particular models considered are, in principle, ruledout by local constraints.Sec. III A begins with reviewing the linear growth ofstructure in Jordan-Brans-Dicke gravity and the result-ing linear matter power spectrum. The description ofcosmological structure in the nonlinear regime is thendiscussed in Sec. III B using N -body simulations and inSec. III C using the spherical collapse model. Sec. III Dreviews different approaches of modelling the chameleoneffect on halo properties such as in the halo mass func-tion and linear halo bias as well as in the density, scalarfield, and mass profiles of the clusters. Finally, Sec. III Efocuses on the description of the nonlinear matter powerspectrum, comparing different modelling approaches, in-cluding nonlinear parametrisations and fitting functions,the halo model decomposition, and perturbation theory. A. Linear perturbations
In scalar-tensor theories, the linear growth functionfor matter fluctuations D ϕ ( a, k ) becomes scale dependentand in the quasistatic limit can be determined from solv- ing D ′′ ϕ + (cid:20) −
32 Ω m ( a ) (cid:21) D ′ ϕ − µ ( a, k )Ω m ( a ) D ϕ ≃ , (23)where Ω m ( a ) ≡ H Ω m a − /H . Hereby, µ ( a, k ) describesthe gravitational modification introduced in the Pois-son equation by the scalar field, which for the Jordan-Brans-Dicke models defined by Eqs. (1) and (2) takesthe form [118, 119] µ ( a, k ) = 1¯ ϕ (cid:20) ω k ¯ ϕa ¯ m + k ¯ ϕ (cid:21) , (24)where ¯ ϕ ≃ m ≃ − α ω (1 − ¯ ϕ ) α − (1 − ¯ ϕ ) α − ¯ R , (25)or ¯ m ≃ [3 f RR ( R = ¯ R )] − in f ( R ) gravity. Eq. (23)can be solved setting the initial conditions D ϕ ( a i , k ) = D ( a i ) = a i and D ′ ϕ (ln a i , k ) = D ′ (ln a i ) = a i at an initialscale factor a i ≪ P L ϕ relates to the ΛCDM power P LΛCDM as P L ϕ ( a, k ) = (cid:18) D ϕ ( a, k ) D ( a ) (cid:19) P LΛCDM ( a, k ) . (26)The ΛCDM variance is S ( a, r ) ≡ σ ( a, r ) = D ( a ) D ( a i ) Z d k | ˜ W ( k r ) | P i ( a i , k ) , (27)where ˜ W ( k r ) is the Fourier transform of a top-hat func-tion of radius r and P i denotes the power spectrum atinitial time a i . Rather than defining the variance at ra-dius r , it can also be written as a function of the mass M = 4 π ¯ ρ m r / λ C ≡ ¯ m − from Eq. (25), Eq. (24) be-comes µ ( a, k ) → D ( a ) in ΛCDM. Note thatwhereas for the chameleon models, Eqs. (23) and (24)strictly only apply in the quasistatic limit, in ΛCDM with µ ( a, k ) = 1, Eq. (23) becomes exact on all scales and canbe derived from combining the linearly perturbed Ein-stein field equations with the energy-momentum conser-vation in the total matter gauge. For the scalar-tensortheories with ΛCDM expansion history considered here,however, deviations in D ϕ ( a, k ) from solving the growthfunction using the full linear perturbation theory aresmall [32, 58, 120–122]. Moreover, the following dis-cussion concentrates on the high-curvature regime whereEqs. (23) and (24) are accurate but it is worth keep-ing in mind that in general, in modified gravity models,corrections to the quasistatic approximation may havemeasurable signatures on horizon scales [58, 123]. B. N -body simulations of chameleon f ( R ) gravity N -body simulations are essential in the study of thechameleon mechanism and its effects on the nonlinearcosmological structure. For this purpose, N -body codeshave been developed by Refs. [85–92]. Thereby, thechameleon suppression has been particularly well stud-ied in the case of the Hu-Sawicki f ( R ) model, for whichin the quasistatic limit, the scalar field and Poisson equa-tions become ∇ δf R = a δR ( f R ) − π G δρ m ] , (28) ∇ Ψ = 16 π G a δρ m − a δR ( f R ) , (29)respectively, where coordinates are comoving, Ψ = δg / (2 g ) is the Newtonian potential in the longitu-dinal gauge, δf R = f R ( R ) − f R ( ¯ R ), δR = R − ¯ R , and δρ m = ρ m − ¯ ρ m . The first simulations solving Eqs. (28)and (29) were performed by Oyaziu et al . [85, 86] witha particle-mesh code on regular grids. The resolution ofthe chameleon force enhancement was then improved us-ing a self-adaptive grid structure by Zhao et al . [89] with mlapm [124] and by Li et al . [125] with ecosmog [90]that is based on ramses [126]. Puchwein et al . [91] intro-duced mggadget , a modification of gadget [127], usingits tree structure to improve the resolution. The equa-tions of motion, Eqs. (28) and (29), in all of these simula-tions are solved in the quasistatic limit, which is a goodapproximation for f ( R ) gravity [85, 92, 122]. Llinaresand Mota [128] introduced a non-static solver for the equations of motion in the symmetron model, which theyapply in its static version in the ramses based N -bodycode isis [92] to simulate the Hu-Sawicki f ( R ) model.The cosmological structure formed in these simula-tions has been analysed through measurements of thematter power spectrum [86, 89–92, 100, 125], halo massfunction [38, 89, 100], halo bias [100], halo densityprofile [46, 99, 100], halo concentration [99], velocitydispersion [47, 74, 97, 99], velocity divergence powerspectrum [125, 129], matter velocity dispersion crosspower [125, 129], redshift-space distortions [125], gravita-tional redshift profiles [130], abundance of subhalos [74],and the integrated Sachs-Wolfe effect [72]. Moreover,hydrodynamical simulations of f ( R ) gravity have beenperformed with mggadget to analyse degeneracies ap-pearing between the signatures of f ( R ) modifications ofgravity and baryonic processes on the matter power spec-trum, also examining the impact from active galactic nu-clei feedback [91]. Further degeneracies with signaturesof massive neutrinos on the simulated matter power spec-trum, halo mass function, and halo bias have been dis-cussed in Ref. [73]. Other studies analysed the hydrody-namics of the intracluster and intragroup medium [74] orcombined simulations with a semi-analytic model for thestudy of galaxy evolution [131].The impact of f ( R ) modifications of gravity on thehalo mass function and matter power spectrum is dis-cussed in Secs. III D and III E. For the numerical com-putations, also involving the spherical collapse density inSec. III C, the cosmological parameters are set to Ω m =1 − Ω Λ with Ω Λ = 0 .
76, Ω b = 0 . h = 0 .
73, slope of the primordial powerspectrum n s = 0 . A s such that the power spectrum normalisationis σ ≡ σ ( a = 1 , r = 8 h − Mpc) = 0 . N -body simu-lations of Ref. [89, 125], from which results for the halomass function and matter power spectrum are also shownin Secs. III D and III E, respectively. The same simula-tion output has been used as comparison for the analyt-ical modelling of these observables in Ref. [75, 99, 103].Here, in order to highlight the effects of modified grav-ity and the chameleon mechanism, numerical computa-tions are specialised to | ¯ ϕ − | = | f R | = 10 − withexponent α = 1 / n = 1), for which the scalar field issmall enough to demonstrate the nonlinear chameleonsuppression effect yet large enough to yield a signifi-cant modification of the cosmological structure. Notethat analogous to Ref. [89], halos will be defined usingthe ΛCDM virial overdensity ∆ vir ≈ f ( R ) halos. Hence, the associated virial masses M vir = 4 π ¯ ρ m ∆ vir r / f ( R ) gravity [100, 103].Finally, it is important to note that Eqs. (28) and (29)form a highly nonlinear system of differential equations,which is harder and takes longer to solve than performing N -body simulations in Newtonian gravity. Additionally,the Hu-Sawicki f ( R ) model is only a subset of scalar-tensor theories or modifications of gravity of cosmologi-cal interest. Hence, it is important to develop simpler,(semi)analytic tools that allow a more general analysis ofthe effects on the nonlinear structure caused by modify-ing gravity. These are ideally also efficient enough for theapplication in Markov Chain Monte Carlo explorations ofthe associated parameter spaces, comparing these com-putations to observations. C. Spherical collapse model
In addition to running N -body simulations ofchameleon gravity, the nonlinear structure formation canalso be studied in the spherical collapse model. Hereby,a dark matter halo is approximated by a constant spher-ically symmetric top-hat density ρ in of radius r th , that isembedded in an outer matter density ρ out . The densitiesare then evolved according to the nonlinear continuityand Euler equations from an initial time to the time ofcollapse of ρ in .In the inner part of the top hat and outside of it, thescalar field is in the minimum of the effective Einsteinframe scalar field potential V eff ( φ ) in Eq. (10) with theJordan frame scalar field values ϕ in and ϕ out denotingthe corresponding solutions for U ϕ ≃ R/
2, for instance,given by Eq. (12). Khoury and Weltman [3] estimate thedistance ∆ r ≥ ϕ ≃ ϕ out to ϕ in as∆ rr th ≃ (3 + 2 ω ) ϕ in − ϕ out κ ρ in r . (30)The force enhancement ∆ F ≡ F − F N at r th due to thenonminimally coupled scalar field, relative to the corre-sponding Newtonian force F N = G m t M/r with testmass m t , can then be approximated by [75, 102, 103]∆ FF N ≃
13 + 2 ω min " rr th − (cid:18) ∆ rr th (cid:19) + (cid:18) ∆ rr th (cid:19) , . (31)This result follows from computing the intermediatescalar field within r ∈ [ r , r th ] in the thin-shell regime∆ r = r th − r ≪ r th , which reproduces the force en-hancement in the thick-shell regime ∆ r > r th when r ≪ r th . Hence, Eq. (31) gives an interpolation be-tween the regime of chameleon suppression ∆ F = 0 andthe regime where the total force F is maximally modi-fied, ∆ F/F N = (3 + 2 ω ) − , which is C for ∆ r/r th → C for ∆ r/r th →
1. Note that the force modificationin Eq. (31) is both dependent on the mass of the halo, M ≡ π ρ in r /
3, and its environmental density ρ out .Li and Efstathiou [102] introduce this force enhance-ment in the evolution of the spherical shell to studythe spherical collapse in chameleon models (cf. [100, 101, 132]). The equation of motion for the physical ra-dius of the spherical top-hat overdensity ζ ( a ) then be-comes [100, 102, 103]¨ ζζ ≃ − κ ρ m − ρ Λ ) − κ (cid:18) FF N (cid:19) δρ m (32)with dots denoting cosmic time derivatives. At initialtime, a i ≪
1, the physical radius of the spherical shellis given by ζ ( a i ) = a i r th but at later times, a > a i , thenonlinear evolution causes ζ ( a ) to deviate from this sim-ple linear relation. This deviation can be characterisedby the dimensionless variable y ≡ ζ ( a ) / ( a r th ), where ρ m / ¯ ρ m = y − follows from conservation of the massenclosed in the overdensity, ¯ ρ m a r = ρ m ζ . Hence,Eq. (32) can be rewritten as y ′′ h + (cid:20) −
32 Ω m ( a ) (cid:21) y ′ h + 12 Ω m ( a ) (cid:18) FF N (cid:19) (cid:0) y − − (cid:1) y h = 0(33)with ρ in / ¯ ρ m = y − and y ′′ env + (cid:20) −
32 Ω m ( a ) (cid:21) y ′ env + 12 Ω m ( a ) (cid:0) y − − (cid:1) y env = 0 , (34)for the environment with ρ out / ¯ ρ m = y − , which follows aΛCDM evolution, corresponding to Eq. (32) in the limitof ∆ F →
0. The force enhancement in Eq. (33) is ob-tained from Eq. (31), replacing ∆ r/r th → ∆ ζ/ζ , wherefor the chameleon models defined by Eq. (13),∆ ζζ ≃ (3 + 2 ω )( ¯ ϕ − a − α − α m ( H r th ) y h Ω Λ Ω m y − + 4 Ω Λ Ω m a ! − α − Ω Λ Ω m y − + 4 Ω Λ Ω m a ! − α . (35)This system of differential equations is then solved withthe initial conditions y h / env , i = 1 − δ h / env , i , y ′ h / env , i = − δ h / env , i , (36)set at an initial scale factor a i ≪ D ϕ ( a, k ) (seeSec. III A), the initial overdensities δ h , i and δ env , i can beextrapolated to a > a i using the ΛCDM linear growthfunction D ( a ) in Sec. III A to define an effective linearoverdensity δ h / env ( x ; ζ h / env ) ≡ D ( a ) D ( a i ) δ h / env , i . (37)In specific, using Eq. (37) to extrapolate the initial over-densities yielding collapse at a given redshift, defines theeffective linear collapse density δ c and the environmen-tal density δ env . Analogous to these two definitions, in @ M Ÿ (cid:144) h D ∆ c X ∆ c H M, ∆ env L\ env ∆ i H r L profiles Σ H M L PPF fitHu - Sawicki È f R0 È = - Α = @ M Ÿ (cid:144) h D ∆ c Designer modelHu - Sawicki model
Α =
Α =
Α =
Α = È f R0 È = - FIG. 2.
Left panel:
Comparison between different approaches to compute the spherical collapse density at z = 0 in theHu-Sawicki f ( R ) model with α = 0 . n = 1) and | f R | = 10 − : environmental average of collapse densities from the thin-shellforce implementation in the spherical collapse in Eqs. (31) and (35) (solid curve); evolution of initial overdensity profile (dashedcurve); reconstructed collapse density from a phenomenological chameleon transition in the variance (dot-dashed curve). Notethat in the mass and environment dependent spherical collapse, the coefficient in Eq. (35) depends on | f R | /M / , hence, δ c canbe scaled to other values of | f R | by a redefinition of mass. Right panel:
Comparison between the spherical collapse densities h δ c i env at z = 0 in the designer (solid curve) and Hu-Sawicki (dashed curves) f ( R ) model for different values α . Eq. (27), the variance S has been defined by the inte-gration of the initial matter power spectrum extrapo-lated to a using the ΛCDM linear growth function D ( a ).Thus, the peak threshold, which will be of interest inSec. III D 1, is determined by ν = δ c / √ S = δ c, i / √ S i = ν i due to the scale-independent growth of structure inΛCDM. In contrast, if using the scale-dependent lineargrowth function D ϕ ( a, k ) of scalar-tensor theories dis-cussed in Sec. III A for these extrapolations instead, thethresholds differ in general, ν = ν i . Hence, Eq. (37) cor-responds to defining the peak-threshold at a i .The gravitational force enhancement determined byEqs. (31) and (35) and thus, δ c , depend on the envi-ronmental density δ env or δ env , i with larger modificationsand stronger suppression for low and high values of δ env ,respectively. In Ref. [102] this environment was specifiedby its Lagrangian (or initial comoving) size with the en-vironment defined as a spherical region around the samecentre as the top-hat overdensity that has a radius thatis larger than the halo that will form but small enough tobe considered its surrounding. Refs. [104, 133] then char-acterised the environment in terms of the Eulerian (phys-ical) size, emphasising that the difference in the respec-tive probability distributions of δ env leads to differencesin the corresponding structures that are formed with theEulerian environments being more likely to be largerthan their Lagrangian counterparts and hence causinga stronger suppression of the modified gravity effects.They argued that a scale at the order of the Comptonwavelength λ C is a natural choice for the Eulerian ra- dius of the environment, setting ζ = 5 h − Mpc. Thisvalue has also been adopted by Refs. [75, 103] togetherwith the probability distribution P ζ ( δ env ) of the Eulerianenvironmental density δ env described in Refs. [104, 134]to analyse the environmental dependence of structuresformed in chameleon theories, also employing differentaveraging procedures. Thereby, using the environmen-tal average of collapse densities, h δ c i env , when modellinghalo properties (see Sec. III D) produced a good match toresults from f ( R ) N -body simulations, furthermore, pro-viding a simplification over modelling observables with δ c ( δ env ) first with subsequent averaging over P ζ ( δ env ).This approach can be further simplified by using the en-vironmental density for which δ c comes close to this av-erage. For the cosmology in Sec. III B and setting z f = 0,one can adopt h δ c i env ≈ δ c ( δ env = 0 . δ env ≈ . h δ env i env ≈ .
16 for the same cosmological configurationhave also been explored in Refs. [75, 103]. The collapsedensity δ c ( δ env = 0 .
4) for collapse at z = 0 and the cos-mology defined in Sec. III B is shown in Fig. 2.Note that there are alternative approaches to comput-ing δ c through the thin-shell force enhancement in thecollapse as described in Eqs. (31) and (35). Schmidt et al . [100] studied the spherical collapse in the extremecases ∆ F/F N = 1 / F/F N = 0 in Eq. (32), cor-responding to the full modification and fully screenedcase in f ( R ) gravity. Borisov et al . [101] generalised thisapproach by considering the isotropic evolution of an ini-tial overdensity profile according to Eqs. (28) and (29),which in this case becomes a one-dimensional problem.They find that, whereas in the limiting cases studied byRef. [100], an initial top-hat overdensity remains a tophat at later times, under this evolution, the initial top-hat overdensity develops a spike at its edge, indicatingshell crossing. To evade the consequent numerical prob-lems, they introduce a Gaussian smoothing of the initialprofile to its cosmological background with its dispersionas a free parameter. However, due to the breakdown ofBirkhoff’s theorem in f ( R ) gravity and the associatedenvironmental dependence, the collapse density becomesdependent on the choice of initial profile and, hence, thedispersion parameter. Kopp et al . [106] followed this ap-proach, introducing a method based on peaks theory andthe matter transfer function that, given the cosmologicalparameters, fixes the initial overdensity profile. More-over, whereas Ref. [101] use the scale and time dependentlinear growth function of f ( R ) gravity in Fourier space,convolved with the Fourier image of a top hat, to define δ c , Ref. [106] define the effective δ c evolved from the ini-tial overdensity according to the linear ΛCDM growthfunction D ( a ) as in Eq. (37). Ref. [106] give a fittingfunction for their δ c as a function of redshift, halo mass,Ω m , and | f R | . The corresponding limit for the sphericalcollapse density in the cosmology specified in Sec. III Bat z f = 0 is shown in Fig. 2. Finally, the chameleon tran-sition in the peak threshold ν can also be modelled viathe phenomenological approach of Li and Hu [95], whichwill be discussed in more detail in Sec. III D 1, wherethe effective collapse density can be reconstructed from δ c ≡ δ Λc p S/S
PPF with δ Λc being the ΛCDM sphericalcollapse density and S PPF is given by Eq. (40).
D. Chameleon clusters
Galaxy clusters are of particular interest when search-ing for cosmological signatures of a chameleon field andplacing constraints on the models. The enhanced growthof structure due to the extra force exerted by the scalarfield yields, for instance, an increase in the abundance ofmassive clusters, as discussed in Sec. III D 1. Counteract-ing this modification is a decrease of the growth enhance-ment near the Compton wavelength of the backgroundscalar field, beyond which gravity returns to Newtonian,causing a decrease of the overabundance of massive ha-los for decreasing values of 1 − ¯ ϕ . More importantly,however, the chameleon mechanism yields a more effec-tive, nonlinear recovery of the Newtonian results, whichis not just dependent on the mass of the halo but alsoon its environment. Similarly, with the increased abun-dance of massive halos, they become less biased, wherethe chameleon mechanism acts again to recover the New-tonian results.Equally interesting are the chameleon field and mat-ter density profiles within the cluster. The gravitationalmodifications yield an increase in halo concentrations for halo masses defined at equal overdensities, accordingly,with enhanced characteristic matter densities and fur-thermore, a matter pile-up in the infall region of the clus-ters. While for values of ϕ ≃
1, light deflection by thismatter distribution is equivalent to the general relativis-tic effect, the dynamically inferred matter distributionsdiffer for unshielded clusters and can be determined fromthe corresponding chameleon field profile as discussed inSecs. III D 2, III D 3, and III D 4.
1. Halo mass function and linear halo bias
The statistics of virialised clusters can be describedusing excursion set theory, where collapsed structures areassociated with regions where the initial matter densityfields smoothed over these regions exceed the threshold δ c . The variance S in Eq. (27) characterises the size ofsuch a region and when varied causes incremental stepsin the smoothed initial overdensity field. These steps areindependent of their previous values if the wavenumbersare uncorrelated, describing a Brownian motion of thesmoothed matter density field with Gaussian probabilitydistribution and with the increment being a Gaussianfield with zero mean. The distribution f of the Brownianmotion trajectories which cross the flat barrier δ c the firsttime at S was characterised by Press and Schechter [135].In the case of chameleon gravity, the barrier is, how-ever, no longer flat and depends on both the variance andthe environment that embeds the collapsing halo. Evenin the Newtonian case, the barrier becomes dependenton the variance once relaxing the assumption of spheric-ity of the halo. Sheth and Tormen [136] modified thePress-Schechter first-crossing distribution based on ex-cursion set results with the moving barrier of ellipsoidalcollapse [137, 138], ν f ( ν ) = N r π q ν h (cid:0) q ν (cid:1) − p i e − q ν / (38)with peak-threshold ν ≡ δ c / √ S , normalisation N suchthat R d ν f ( ν ) = 1, as well as p = 0 . q = 0 . q is set to match the halo mass function n ln M vir ≡ d n d ln M vir = ¯ ρ m M vir f ( ν ) d ν d ln M vir (39)with measurements from ΛCDM N -body simulations.Although Eq. (38) has been obtained in the concordancemodel, in combination with the mass and environmentdependent spherical collapse model in Sec. III C, it alsoprovides a good description of the relative enhancementof the chameleon halo mass function with respect toits ΛCDM counterpart using the same values of p and q [75, 103]. The corresponding enhancements for the Hu-Sawicki and designer f ( R ) models are shown in Fig. 3.Also shown are the enhancements measured in the Hu-Sawicki α = 0 . N -body simulations of Ref. [125] (black0 @ M Ÿ (cid:144) h D D n (cid:144) n X ∆ c H M, ∆ env L\ env ∆ i H r L profiles Σ H M L PPF fitHu - Sawicki È f R0 È = - Α = @ M Ÿ (cid:144) h D D n (cid:144) n Designer modelHu - Sawicki model
Α =
Α =
Α =
Α = È f R0 È = - FIG. 3.
Left panel:
Relative enhancement in the halo mass function at z = 0 of the Hu-Sawicki | f R | = 10 − model ( α = 0 . N -body simulations of Ref. [125] (black data points) and Ref. [89] (gray data points). Rightpanel:
Comparison between the corresponding enhancements in the designer model and the Hu-Sawicki model with differentvalues of α . data points) and Ref. [89] (gray data points). In com-parison, the designer model yields a smaller increase inthe halo mass function for massive halos, whereas theenhancement is larger for small masses. This behaviourreflects the differences in δ c and f R /f R between the twomodels illustrated in Figs. 1 and 2.Instead of using the phenomenological halo mass func-tion in Eq. (38), one can use excursion set theory to com-pute a theoretically better motivated first-crossing dis-tribution with the moving barrier defined by the linearchameleon collapse density δ c ( S, δ env ). The correspond-ing halo mass function is then determined from Eq. (39).This approach was conducted in Ref. [102] using a La-grangian definition of environment and extended to Eu-lerian environments in Ref. [104], who compared the twoapproaches as well as performed both numerical integra-tions and Monte Carlo simulations, finding that due tothe larger likelihood of high-density environments in theEulerian case, the overabundance of medium and largesize halos with respect to ΛCDM is weakened. The halomass function of the Hu-Sawicki f ( R ) model using ex-cursion set theory and a numerical integration methodhas been computed in Ref. [103]. The comparison to N -body simulations showed a better agreement with theSheth-Tormen approach when combined with the massand environment dependent spherical collapse model anda subsequent averaging over the probability distributionof the Eulerian environment. Ref. [106] formulated ananalytic expression for the halo mass function based onexcursion set theory with a drifting and diffusing barriercomputed from the f ( R ) evolution of the initial densityprofile with uncorrelated steps, which they tested against Monte Carlo random walks. Note that in Fig. 3, the re-sults referring to the smoothed initial density profile ofRef. [106] are not using this result for the halo mass func-tion but apply the spherical collapse density from thisapproach (see Sec. III C) to the Sheth-Tormen formulain Eq. (38).Finally, whereas in all of these approaches, the ran-dom walk was considered Markovian, i.e., with uncorre-lated steps, Lam and Li [105] introduce correlated stepsin the excursion set approach and find that, in general,this leads to an enhancement of the modifications in thehalo mass function. The difference with respect to uncor-related steps is due to a change of distribution of the Eu-lerian environmental densities, whereas the Lagrangiandefinition is not affected, and hence, the first-crossingprobability, as well as correlations between δ and δ env .In order to analyse these effects in the first-crossing dis-tribution, they study three different window functions, asharp k -filter for the uncorrelated case, as well as a Gaus-sian and top-hat filter for the correlated steps, for whichthey run Monte Carlo simulations.Alternatively to introducing the chameleon mechanismin the spherical collapse density, Li and Hu [95] proposeda nonlinear parametrised post-Friedmann (PPF) descrip-tion to determine the halo mass function for chameleon f ( R ) gravity. They phenomenologically interpolate be-tween the linearised and suppressed regimes by introduc-ing a chameleon PPF transition in the variance as S / ( M ) = S / ϕ ( M ) + ( M/M th ) µ S / ( M )1 + ( M/M th ) µ , (40)where M th and µ are fitted simultaneously to different1halo mass functions extracted from N -body simulationswith different configurations of the chameleon field. ThePPF peak threshold in Ref. [95] is then given by ν PPF ≡ δ Λc S / ( M ) , (41)which they subsequently use in the Sheth-Tormen ex-pression Eq. (38) to approximate the halo mass func-tion. The parameters ( M th , µ ) have been calibratedto N -body simulations of the Hu-Sawicki f ( R ) grav-ity model in Refs. [95, 103]. With the mass definition M from setting ∆ = 200, Ref. [95] finds M th =1 . × (cid:0) | f R | (cid:1) / M ⊙ /h and µ = 2 . M following Ref. [103]. Note that if attributing thechameleon transition in the peak threshold ν PPF and ac-cordingly in the halo mass function to a modification ofthe spherical collapse density rather than to a transitionin the variance (see Sec. III C), the scaling of M th with | f R | / can be derived from the spherical collapse modelwhen considering the coefficient of the force modificationin Eq. (35) that scales as (1 − ϕ ) /r ∼ (1 − ϕ ) /M / .Finally, with the Sheth-Tormen halo mass functionEq. (38), the linear halo bias obtained in the peak-background split becomes [136] b L ( M vir ) ≡ b ( k = 0 , M vir ) = 1+ a ν − δ c + 2 pδ c [1 + ( a ν ) p ] . (42)The effective linear collapse density δ c in chameleon mod-els is suppressed relative to ΛCDM, which causes b L todecrease. Moreover, the modification becomes mass andenvironment dependent and can be determined using thespherical collapse model described in Sec. III C.Ref. [100] performed a measurement of the linear halobias in N -body simulations of Hu-Sawicki f ( R ) grav-ity and found good agreement with the deviations in b L predicted by Eq. (42) when using the modified peak-threshold. Halo biasing has also been analysed forchameleon models in Ref. [105] within the framework ofexcursion set theory using the unconditional and condi-tional first-crossing distribution with different smoothingwindow functions. They find an increase in the modifica-tions of b L for correlated steps with respect to the uncor-related case due to less likely high-density environmentsand correlations of δ env with δ in the non-Markovian sce-nario. Deviations in the halo bias and halo mass functionhave also been analysed for a Yukawa-type modificationof gravity, using the spherical collapse model and excur-sion set theory based on the scale-dependent modificationof the growth function [139–141].
2. Cluster density profile
Navarro, Frenk, and White (NFW) [142] found thatthe dark matter clusters formed in ΛCDM N -body sim- ulations are well described by spherical halos with thesimple universal radial density profile δρ m ( r ) = ρ s rr s (cid:16) rr s (cid:17) . (43)Hereby, ρ s and r s denote the characteristic density andscale, respectively, which can be calibrated to simula-tions. Alternatively, given a specific virial halo mass M vir that is defined by the virial overdensity ∆ vir , one may usethe virial halo concentration c vir ≡ r vir /r s as the free pa-rameter and the relations ρ s = 13 ¯ ρ m ∆ vir c (cid:20) ln(1 + c vir ) − c vir c vir + 1 (cid:21) − , (44) r s = 1 c vir (cid:18) M vir π ¯ ρ m ∆ vir (cid:19) / , (45)in the density profile Eq. (43). The NFW profile canfurther be reduced to a function of halo mass only byadopting a mass-concentration scaling relation. In thefollowing, the concentration shall be given by the relation c vir ( M vir , a ) = 9 a (cid:18) M vir M ∗ (cid:19) − . (46)with the critical mass M ∗ satisfying S ( M ∗ ) = δ .Eq. (46) has been calibrated to approximately 5 × halos of mass 10 − M ⊙ /h extracted from ΛCDM N -body simulations in Ref. [143] and shall here be as-sumed to also apply to more massive halos, concentratingon ∆ vir ≡
390 as in Sec. III B.When gravity is modified, Eqs. (43) and (46) maynot be applicable and it is worth checking their perfor-mance against N -body simulations. Ref. [99] found thatEq. (43) provides comparably good fits to the dark mat-ter halo density profiles extracted from N -body simula-tions of the Hu-Sawicki f ( R ) model as to the concor-dance model halos. Given that the fit is accurate fora range of | f R | values and, accordingly, for differentmagnitudes of the gravitational modifications and thesimilarity of the chameleon mechanism between differ-ent choices of model parameters in Eq. (35), this moti-vates the use of the NFW profiles for a wider class ofchameleon models defined by Eqs. (1), (2), and (13).The scaling relation Eq. (46) can be adopted in thechameleon models with the replacement S → S ϕ , where S ϕ is the variance computed using the scale-dependentgrowth function D ϕ ( a, k ) [46, 95, 100]. In this approach,however, the chameleon modification is incorporated in M ∗ in a linear manner, which depends on the ampli-tude ( ¯ ϕ −
1) and only rescales c vir ( M vir ) with a con-stant factor. Hence, it does not capture the chameleonmechanism, neither taking into account dependencies ofthe modification on mass nor environment [99]. Thesedependencies can be introduced in c vir by reinterpret-ing the concentration-mass relation defining the criticalmass as M ∗ ( δ c , σ ) ≡ σ − ◦ δ c and adopting the effective2 δ c from the chameleon spherical collapse model [75]. Thechameleon mechanism is then incorporated in the con-centration via the relations c vir ( M vir , δ env , a ) = 9 a (cid:20) M ∗ ( M vir , δ env ) M vir (cid:21) . , (47) M ∗ ( M vir , δ env ) ≡ ( σ − ◦ δ c )( M vir , δ env )= σ − ( δ c ( M vir , δ env )) . (48)As δ c becomes smaller in the chameleon models, M ∗ andthe concentration become larger, causing an enhance-ment in ρ s and a decrease of r s compared to their ΛCDMcounterparts. This behaviour is in agreement with themeasurements of halo concentration, characteristic den-sity, and characteristic radius from N -body simulationsof the Hu-Sawicki f ( R ) model [99].Finally, while within the virial radius, the halo densityprofiles of the chameleon and concordance model clustersagree up to different values of characteristic density andradius, at a few virial radii, N -body simulations showan enhancement of the halo-matter correlation functiondue to a matter pile-up in the infall region caused bythe late-time enhanced gravitational forces [46, 100] (alsosee [71]). The effect can be well described by the halomodel [46].
3. Chameleon field profile
Given the dark matter profile of a cluster discussed inSec. III D 2, i.e., the NFW relation for δρ m in Eq. (43),the scalar field profile within the halo can be obtainedfrom solving the scalar field equation, Eq. (10). As-suming sphericity and the quasistatic limit, this yieldsa second-order differential equation for ϕ ( r ), which caneasily be integrated numerically adopting the substitu-tion ϕ − − e u ( r ) [85, 89, 97, 99].Alternatively, an analytic approximation for thechameleon field profile can be derived requiring ϕ ≃ U ( ϕ ) in Eq. (10)with respect to ¯ ϕ . After subtraction of the background, δϕ ≡ ϕ − ¯ ϕ , the scalar field equation becomes˜ ∇ δϕ lin − ¯ m δϕ lin + κ ω δρ m ≃ , (49)which is solved by [99] δϕ lin ≃ − κ ρ s r ω n Γ[0 , ¯ m ( r + r s )] e m ( r + r s ) +Γ[0 , − ¯ m ( r + r s )] − Γ(0 , − ¯ m r s ) − e m r s Γ(0 , ¯ m r s ) (cid:9) e − ¯ m ( r + r s ) r (50)with the upper incomplete gamma functionΓ( s, r ) = Z ∞ r d t t s − e − t . (51) The integration constants are set by the requirementsthat lim r →∞ δϕ lin = 0 and lim r → rδϕ lin = 0. In thisapproximation, the chameleon transition is assumed totake place instantaneously once δϕ = 1 − ¯ ϕ , hence, δϕ ≈ min ( δϕ lin , − ¯ ϕ ) (52)or equivalently, ϕ ≈ min ( ϕ lin , ρ m ≫ ¯ ρ m , thelinearised scalar field simplifies to [99] δϕ lin ≃ κ ρ s r ω (cid:20) ln(1 + r/r s ) r − ¯ m e ¯ m r s Γ(0 , ¯ m r s ) (cid:21) , (53)for which the chameleon screening scale is [75] r c = − r s − A − W [ − A r s exp( − A r s )] (54)with the Lambert W function W [ · ], solving x = W ( x ) exp[ W ( x )], and A ≡ ωκ ρ s r (1 − ¯ ϕ ) + ¯ m e ¯ m r s Γ(0 , ¯ m r s ) . (55)Rather than assuming an instantaneous transitionof the linearised scalar field into a chameleon-shieldedregime, one may want to require a continuously differen-tiable transition in δϕ [98]. In this case, a free r c and thetwo integration constants of the outer linearised solution δϕ out obtained from the integration of Eq. (49) in thelimit of ρ m ≫ ¯ ρ m , i.e., neglecting the term − ¯ m δϕ , arematched to the inner chameleon-shielded solution δϕ in ,where U ϕ ≃ R/
2, which for the potential Eq. (13) takesthe expression Eq. (12). More specifically, the chameleonfield δϕ = (cid:26) δϕ out , r > r c δϕ in , r ≤ r c (56)and r c are required to satisfy the conditions δϕ out ( r c ) = δϕ in ( r c ) , (57) δϕ ′ out ( r c ) = δϕ ′ in ( r c ) , (58)and that lim r →∞ ϕ out ( r ) = ϕ env , where ϕ env ≃ ¯ ϕ willbe assumed. While in the resulting relations, one inte-gration constant is set by the environment, the other isset by the transition scale r c , which has to be computednumerically.The different profiles discussed here have been com-pared to each other in Ref. [99] showing good agree-ment with N -body simulations of the Hu-Sawicki f ( R )model. Note, however, that in all of these approaches,requiring that lim r →∞ δϕ out ( r ) = 0 does not necessarilyyield a recovery of the simulated chameleon field profileof a particular cluster as the environment of the clus-ter can deviate from the cosmological background. Inthis case, δϕ should be matched to boundary conditionssuch as obtained from extracting δϕ ( r vir ) from simula-tions [99]. Furthermore, note that the scalar field profileobtained from Eqs. (56), (57), and (58) does not apply3on scales approaching the Compton wavelength, wherefor scales beyond it the scalar field is suppressed accord-ing to Eq. (50).Following Refs. [48, 77], instead of solving r c numeri-cally in the last approach, given that ρ m ≫ ¯ ρ m , one canapproximate the inner solution in Eq. (56) as ϕ in ≃ r c ≃ κ ρ s r ω − ¯ ϕ − r s , (59) δϕ out ≃ κ ρ s r ω ln (cid:18) r + r s r c + r s (cid:19) r + (1 − ¯ ϕ ) r c r . (60)
4. Dynamical mass profile
An unscreened test particle of mass m t in the clusterexperiences the extra force F ϕ ≡ − m t ∇ ϕ. (61)due to the presence of the chameleon field. This extraforce can be interpreted as being exerted by the phantommass M ϕ ≡ − r G dd r ϕ ( r ) , (62)assuming a spherical system and ϕ ≃
1. Eq. (62)contributes to the dynamically inferred mass as M D = M + M ϕ . Using Eqs. (56), (59), and (60), it thereforefollows that M D ( r ) ≃ (cid:26) r − r c )3 + 2 ω (cid:20) − M ( r c ) M ( r ) (cid:21)(cid:27) M ( r ) , (63)where Θ is the Heaviside step function. The cluster massprofile M ( r ) = 4 πρ s r (cid:20) ln (cid:18) rr s (cid:19) − rr + r s (cid:21) + M c (64)is obtained from integrating the NFW density profile,where M c is a mass correction due to deviations fromthis relation in the inner part of the cluster. The approx-imations and equations leading to the expression Eq. (63)generalise and yield a derivation of the assumption madefor M D ( r ) in Ref. [97]. Note, however, that Eq. (63)does not apply on scales approaching the Compton wave-length, where for scales beyond it the scalar field is sup-pressed according to Eq. (50) (see Secs. III A and III D 3).The difference between M D ( r ) and M ( r ) also depends onthe environment, which in Eq. (63) can be taken into ac-count through dropping the condition ϕ env = ¯ ϕ . It hasalso been shown in Ref. [133] that at r = r vir , the dif-ference in mass observed in N -body simulations of f ( R )gravity can be well described using the thin-shell condi-tion.Moreover, dense objects such as stars may not feel thefull force modification in Eq. (63) due to self-shielding in the chameleon mechanism. The effective force felt bysuch an object can be estimated by Eqs. (30) and (31),where ϕ out is determined by its environmental field value.While the stars may be screened or partially screened, thegas of a cluster feels the full force modification. Assum-ing hydrostatic equilibrium of the gas, for a sphericallysymmetric system, the gas density ρ gas and pressure P relate to the dynamical mass profile as1 ρ gas ( r ) d P ( r )d r = − G M D ( r ) r . (65)Assuming no contribution from non-thermal pressure,the gas pressure relates to the gas density and tempera-ture as P = P thermal ∝ ρ gas T gas . Whereas the mass pro-file M ( r ) = M D ( r ) − M ϕ ( r ) can be measured by weakgravitational lensing around a cluster, the gas tempera-ture, density, and pressure can be determined from X-ray observations and the Sunyaev-Zel’dovich effect. Inhydrostatic equilibrium, these observables are uniquelydetermined from any combination of two of these pro-files. Hence, a combination of these measurements yieldsa powerful test of gravity [77]. E. Matter power spectrum
In the halo model [144–146], statistics of cosmologi-cal structures are decomposed into the underlying halocontributions. The nonlinear matter power spectrum isdescribed by the two-halo and one-halo terms, P mm ( k ) ≃ P h ( k ) + P h ( k ) , (66) P h ( k ) = Z d ln M vir n ln M vir M ¯ ρ | y ( k, M vir ) | , (67)where y ( k, M ) denotes the Fourier transform of the halodensity profile, which shall be given by the NFW expres-sion Eq. (43), with a truncation at r vir and normalisationlim k → y ( k, M ) = 1. The two-halo term shall be approx-imated by P h ( k ) ≈ P L ( k ), and the halo mass functionand halo concentration shall be determined as describedin Secs. III D 1 and III D 2. Note that the use of the lin-ear matter power spectrum as the two-halo contributionunderestimates nonlinear effects on the transition scalesbetween the one-halo and two-halo terms. Although thisis also a problem for modelling the ΛCDM power spec-trum, for chameleon models, the situation is further com-plicated as the suppression mechanism is not incorpo-rated in the linearly computed growth enhancement. InRef. [75], it was shown that the relative enhancement inthe chameleon matter power spectrum with respect to itsΛCDM counterpart can be well described by introducinga simple scale-dependent correction to the linear powerspectrum P effL ϕ ( a, k ) = P L ϕ ( a, k ) + ( k/k ∗ ) P LΛCDM ( a, k )1 + k/k ∗ , (68)4where k ∗ ≈ . p (1 − ¯ ϕ ) / − h Mpc − . The scale de-pendence in this correction can be motivated by the re-lation between the top-hat size and scalar field ampli-tude in the coefficient of the force modification used inthe spherical collapse computation in Eq. (35). The re-sulting enhancements of the matter power spectrum forthe Hu-Sawicki and designer f ( R ) models are shown inFig. 4, where initial power spectra are determined us-ing the Eisenstein-Hu transfer function [147, 148]. Thehalo mass function and concentration in the one-haloterm Eq. (67) are computed using the different spheri-cal collapse densities described in Sec. III C. In the caseof the mass and environment dependent spherical col-lapse model, the most probable environment is assumed,which corresponds to the environmental density at thepeak of its probability distribution. This can be inter-preted as an averaging procedure over the unshielded,shielded, and partially shielded forces acting on the darkmatter particles.A similar approach to the correction of the transitionregime in Eq. (68) was conducted in Ref. [95], who replacethe two-halo term and its interpolation to the one-haloterm with the same phenomenology as halofit [149].More specifically, P h → π ∆ ( k ) /k with∆ ( k )∆ ( k ) = (cid:2) ( k ) (cid:3) β n α n ∆ ( k ) exp (cid:18) − y − y (cid:19) , (69)where ∆ ( k ) = (2 π ) − k P L ( k ) and y = k/k σ determinethe transition scale to the one-halo term with σ G ( k − σ ) =1 and σ ( R ) = Z d ln k ∆ ( k ) exp( − k R ) . (70)The transition parameters α n and β n are determinedfrom [149] α n = 1 . . n eff − . n , (71) β n = 0 . . n eff + 0 . n , (72)where n eff ≡ − − d ln σ ( R )d ln R (cid:12)(cid:12)(cid:12)(cid:12) σ G =1 . (73)Predictions of P ( k ) then vary depending on whetherthe modified or ΛCDM linear matter power spectrum isused to determine the transition. The arithmetic meanbetween the two power spectra produced by using ei-ther P LΛCDM or P L ϕ to compute the right-hand side ofEq. (69) and σ G is shown in Fig. 4. Note that in this case,the one-halo term is computed according to Ref. [95], us-ing the halo mass function and halo concentration de-termined from the PPF fit in the variance discussed inSec. III D 1.A different nonlinear PPF formalism, directly mod-elling P ( k ), was proposed by Hu and Sawicki [93].It interpolates between the modified and suppressed regimes of the chameleon matter power spectrum usingthe nonlinear power spectrum in ΛCDM P ΛCDM for theshielded regime and in the unshielded regime, its coun-terpart P non − ΛCDM that extrapolates the modified linearpower spectrum to small scales and does not exhibit achameleon transition. The power spectra are then inter-polated via the relation P ( a, k ) = P non − ΛCDM ( a, k ) + c nl Σ ( a, k ) P ΛCDM ( a, k )1 + c nl Σ ( a, k ) . (74)In Ref. [93], the transition was assumed to scale asΣ ( k ) = (2 π ) − k P L ( k ). Koyama et al . [107] used per-turbation theory and Σ ( k ) = (cid:2) (2 π ) − k P L ( k ) (cid:3) / torecover the simulated power spectra in f ( R ) gravity upto k ∼ . h Mpc − . Zhao et al. [94] then extended thePPF approach for P ( k, a ), fittingΣ ( a, k ) = (cid:20) k π P L ( a, k ) (cid:21) α nl + β nl k γ , (75)to N -body simulations of the Hu-Sawicki model with α = 0 .
5, where c nl , α nl , β nl , and γ are re-calibrated atdifferent redshifts and for different values of | f R | . Thecorresponding matter power spectrum using halofit todetermine P ΛCDM and P non − ΛCDM is shown in Fig. 5. Re-cently, Zhao [96] introduced mghalofit , which directlymodifies halofit by generalising its fitting parameters toaccommodate the Hu-Sawicki f ( R ) model. The enhance-ment in P ( k ) predicted by mghalofit is also shown inFig. 5.Alternatively, the linear, quasilinear, and nonlinearscales of the matter power spectrum can be computedcombining perturbation theory at the one-loop level withthe halo model. An implementation of this concept wasdeveloped by Valageas et al . [150] and applied to differ-ent scalar-tensor theories, including the Hu-Sawicki f ( R )model, by Brax and Valageas [108]. In order to computethe one-loop matter power spectrum, the scalar modifi-cation of the Newtonian potential is expanded to thirdorder in the nonlinear density fluctuation. Subsequently,a Lagrangian-space regularisation is applied to the one-loop expansion and combined with the halo model, usingthe gravitational modification on the high-mass end ofthe halo mass function determined by a modified spheri-cal collapse calculation and the NFW density profile withidentical mass-concentration relation to the concordancemodel, hence, not accounting for modified gravity effects.The results obtained in Ref. [108] for the enhancementof the matter power spectrum in the Hu-Sawicki f ( R )model with α = 0 . f ( R ) gravity, nonlinear correc-tions to the matter power spectrum begin to contribute atslightly larger scales than in ΛCDM. The enhancementsin the power spectrum at z = 0 predicted by linear theorycan be up to 50% and 100% larger at k = 0 . h Mpc − and k = 0 . h Mpc − , respectively, than their nonlinearcounterparts.5 @ h (cid:144) Mpc D D P (cid:144) P X ∆ c H M, ∆ env L\ env ∆ i H r L profiles Σ H M L PPF fitHu - Sawicki È f R0 È = - Α = @ h (cid:144) Mpc D D P (cid:144) P Designer modelHu - Sawicki model
Α =
Α =
Α =
Α = È f R0 È = - FIG. 4.
Left panel:
Relative enhancement of the matter power spectrum at z = 0 in the Hu-Sawicki model for α = 0 . | f R | = 10 − with respect to ΛCDM predicted by halo model approaches with the different spherical collapse densities of Fig. 2.The data points represent results from N -body simulations of Ref. [125]. Right panel:
Comparison between the matter powerspectrum enhancements at z = 0 in the designer model (solid curve) and the Hu-Sawicki model (dashed curves) with differentvalues of α using the halo model and the mass and environment dependent spherical collapse model. Dotted curves indicatelinear results. @ h (cid:144) Mpc D D P (cid:144) P P H k L PPF fitMGHalofitLinear P H k L Hu - Sawicki È f R0 È = - Α = - loop (cid:144) halo model FIG. 5. Same as left panel of Fig. 4 but using mghalofit (solid curve), the PPF fit Eq. (74) for P ( k ) (dashed curve), acombination of one-loop perturbation theory and halo model(dot-dashed curve), and linear perturbations (dotted curve). IV. OBSERVATIONAL CONSTRAINTS
Chameleon models have been constrained using a va-riety of observations from laboratory to cosmologicalscales [8, 50, 59, 65, 75, 77, 82, 83, 113, 151]. Fig. 1shows the current bounds inferred on the field ampli- tude and coupling strength on local [8, 75], astrophysi-cal [50], and cosmological scales [77]. In particular, theHu-Sawicki and designer f ( R ) models ( ω = 0) have beentested using a range of observables and methods. A selec-tion of current and prospective constraints on | f R | aresummarised in Table I, focusing mainly on cosmologicalresults. Note that due to the difference in the signa-tures of the Hu-Sawicki and designer f ( R ) models, forrigour, the 95% confidence level constraints listed fromthe different analyses are assigned to the specific formof f ( R ) assumed, i.e., the designer (D) and Hu-Sawicki(HS) models. Hereby, the Hu-Sawicki case refers to themodel specifications with α = 0 . n = 1). A fewof the results listed in Table I shall shortly be reviewedhere.With increasing | f R | , the growth of structure becomesenhanced at late times. This affects the cosmic mi-crowave background (CMB) through modifications of theintegrated Sachs-Wolfe (ISW) effect, gravitational lens-ing, and the Sunyaev-Zel’dovich (SZ) effect. With grow-ing | f R | , the ISW contribution is reduced, initially yield-ing a decrease of the temperature anisotropy power spec-trum at low multipoles. Eventually, with increasing mod-ification, the ISW contribution to the temperature fieldbecomes negative, at which point, however, the ISW con-tribution to the temperature anisotropy power spectrumstarts to rise again due to the square in the temperaturefield. While the initial reduction of the ISW contributionis slightly preferred by the data, the enhancement with | f R | after the turning point can be used to place con-straints on the model at the order of | f R | . − [23, 32].6 TABLE I. Selection of current and prospective constraints on f ( R ) gravity.Measurement redshift scale | f R | constraint model Ref.Integrated Sachs-Wolfe (ISW) effect . & h − Mpc < . × − D [23, 32]Galaxy-ISW cross correlations . & h − Mpc < . × − D [9, 30, 32]Galaxy power spectrum (WiggleZ) 0 . − & h − Mpc 10 − . | f R | . × − D [79]Galaxy power spectrum (WiggleZ) 0 . − & h − Mpc . . × − D [79]Redshift-space distortions (LRG) 0 . − .
47 (15 − h − Mpc . − HS [33] E G probe 0 .
32 (10 − h − Mpc .
10 D [31, 32]CMB lensing (ACT) . − h − Mpc . − D [60]CMB lensing (SPT) . − h − Mpc . × − D [60]CMB lensing (Planck) . & h − Mpc . − D [66, 68]Cluster abundance (Chandra) < .
15 (1 − h − Mpc < . × − HS [29, 38]Cluster abundance (MaxBCG) 0 . , .
25 (1 − h − Mpc < . × − D [32]Gravitational redshift of galaxies in clusters 0 . − .
55 (0 . − h − Mpc . . . HS [45]Cluster density profiles (maxBCG) 0 .
23 (0 . − h − Mpc < . × − HS [46]Coma gas measurements 0 .
02 (0 . − h − Mpc < × − D/HS [77]Strong gravitational lenses (SLACS) 0 . − .
36 (1 −
10) kpc < . × − HS [28]Solar System 0 .
20 au / < × − D/HS [8, 75]Supernova monopole radiation ∼ ∼ R ⊙ . − D/HS [65]Distance indicators in dwarf galaxies . . . R ⊙ < × − D/HS [50]Relativistic effects in galaxy-clustering . & h − Mpc . − D [58]21 cm intensity mapping + CMB 0 . − . & h − Mpc . − D [55]CMB ISW-lensing bispectrum . & h − Mpc . − D/HS [53, 81]Matter bispectrum ∼ . h − Mpc 10 − . | f R | . − HS [44]Stacked phase-space distribution 0 . − . − h − Mpc . (10 − − − ) HS [47]Galaxy infall kinematics 0 .
25 (0 . − h − Mpc . (10 − − − ) HS [71]Dwarf galaxies ∼ ∼ . − D/HS [41, 61]
When cross correlating the ISW effect with foregroundgalaxies, the linear contribution of the temperature fieldyields an anti-correlation for strong modifications, whichcan be used to tighten constraints on | f R | by about a fac-tor of 5 [30, 32]. With the improved measurement of thehigh angular multipoles, especially by the Planck mis-sion, it has also become possible to constrain | f R | withCMB lensing, tightening constraints by almost one orderof magnitude over the galaxy-ISW bounds [60, 66, 68].Using the galaxy power spectrum measured by Wig-gleZ , Dossett et al . [79] recently derived a constraint of10 − . | f R | . × − and | f R | . − on the de-signer model, assuming that the modifications at z ∼ . k = 0 . h Mpc − and k = 0 . h Mpc − ,respectively. Measurements of the galaxy clustering in asample of luminous red galaxies (LRG) from the SloanDigital Sky Survey (SDSS) combined with the galaxy-galaxy lensing signal and galaxy velocities obtained fromredshift-space distortions were used by Reyes et al . [31]to get a measurement of the E G parameter [22] thatprovides a robust test of gravity and cancels uncertain-ties in galaxy bias and the initial amplitude of matterfluctuations, however, only yielding weak constraints on | f R | [32].The enhanced abundance of massive clusters discussedin Sec. III D 1 has been used by Schmidt et al . [29], using Chandra
X-ray data, and Ref. [32], using SDSS MaxBCGclusters, to constrain the Hu-Sawicki ( α = 0 .
5) and thedesigner model, respectively, at the level of | f R | . − .A prescription for mapping these constraints to differ-ent values of α in the Hu-Sawicki model was formu-lated by Ferraro et al . [38]. Similarly, the enhancedabundance affects the stacking of cluster density profiles,which furthermore, shows a signature of a matter pile-upin the infall region due to the late-time enhanced gravi-tational forces in f ( R ) gravity. Weak gravitational lens-ing measurements around maxBCG clusters were usedin Ref. [99] to constrain these effects and place an up-per bound of | f R | . − on the model. Comparableconstraints to the ones inferred from the increased abun-dance have been estimated by Yamamoto et al . [33] forredshift-space distortions measured in the LRG sample.Interesting constraints on modified gravity can alsobe obtained using the difference in dynamically inferredmasses to masses inferred from weak lensing discussed inSec. III D 4. Refs. [28, 152] constrained modifications ofgravity measuring the dynamical masses of strong lenses7via their stellar velocity dispersions and using the radiiof Einstein rings to determine the lensing masses in asample of elliptical galaxies from the Sloan Lens ACS(SLACS) survey. In specific, Smith [28] derived an up-per bound of | f R | < . × − on the Hu-Sawicki f ( R )model ( α = 0 . f ( R ) gravity using Eq. (63),where massive halos recover the results of Newtoniangravity, and pointed out that this modification is alsoreduced for halos which have more massive halos in theirproximity. This mass estimator was used by Wojtak etal . [45] to constrain modifications of gravity by disentan-gling the gravitational redshift of the light emitted bygalaxies, which propagates through the cluster and shiftsthe observed galaxy centre, from the kinetic Doppler ef-fect due to galaxy motions, which broadens the observedvelocity distribution. Their measurement is compati-ble with the f ( R ) modifications. Recently, Terukina etal . [77] have combined gas measurements in the Comacluster of the X-ray temperature and surface brightnessas well as from the SZ effect with lensing measurementsaround the cluster to infer a constraint on f ( R ) gravityof | f R | . × − .Finally, the currently strongest constraints on f ( R )gravity in astronomy are inferred from local and astro-physical tests of gravity (see Fig. 1). The requirementon the model to satisfy Solar System constraints and theconsequent interpolation of the chameleon field from thishigh-curvature regime to the low curvature on galacticscales and the environment of the Milky Way, assumedto be the cosmological background, yields a constraint of | f R | < × − [8, 75]. Note that for the chameleonmodels given by Eq. (13), local tests of the equivalenceprinciple such as conducted by the lunar laser rangingexperiment, E¨ot-Wash, or the Cassini mission alone onlyplace weak constraints on the modifications of gravity ifassuming that the scalar field in the galactic backgroundis in the minimum of its effective potential. This corre-sponds to assuming that the scalar field in our Galaxyor the Solar System satisfies the high-curvature solutiondiscussed in Sec. II A, in which case the chameleon field issuppressed proportionally to the ratio between the localand background curvature. Constraints from the condi-tion of a locally minimised chameleon field can be in-ferred by requiring that the background scalar field be-comes smaller than the galactic gravitational potentialwithin the characteristic scale of the halo, which can berelated to constraints on maximal rotation velocities [8],or similarly, by modelling the matter distribution of theMilky Way, solving the resulting scalar field equation,and tracing the chameleon field from the cosmologicalbackground to the location of the Solar System, wherethe high-curvature solution requires it to be shielded [75].The comparison of shielded and potentially unshieldeddistance indicators from tip of the red giant branch starsand cepheids, respectively, in a sample of unscreened dwarf galaxies gives a constraint of | f R | < × − [50].Another interesting but weaker astrophysical constraintwith an upper bound of | f R | . − can furthermorebe inferred from the absence of monopole radiation incore-collapse supernovae [65]. V. OUTLOOK
With the increasing wealth and quality of observa-tional data that will be collected with future surveys,constraints on the gravitational models can be improvedand novel methods for testing scalar-tensor gravity andchameleon models will become feasible. Prospective con-straints on f ( R ) gravity from a variety of observables andon a wide range of scales have, for instance, been anal-ysed in Refs. [41, 44, 47, 53, 55, 58, 61, 71, 81]. The pos-sibility of using relativistic corrections in horizon-scalegalaxy clustering, measured in a multi-tracer analysis tocancel cosmic variance, for inferring constraints on darkenergy and modified gravity models, has been exploredin Ref. [58], finding a prospective bound of | f R | . − on the designer f ( R ) model. Hall et al . [55] have esti-mated that the combination of 21 cm mapping and CMBdata will allow to place a constraint on this field ampli-tude of order 10 − and a modification of the same ordershould also be detectable in the bispectrum of the darkmatter density field [44]. The cross correlation of theISW effect with gravitational lensing generates a signa-ture in the bispectrum of the CMB temperature field.With enhanced growth at late-times, f ( R ) gravity sup-presses the ISW-lensing cross correlation and modifiesthe temperature bispectrum. This effect was used by Hu et al . [53] and Munshi et al . [81] to forecast a constraintof | f R | . − for Planck results.The enhanced late-time gravitational forces also in-crease velocities of unshielded objects. Signatures of thiseffect have, for instance, been characterised in the ve-locity power spectrum of dark matter particles [125], forredshift-space distortions [33, 54, 56, 129], in the stackedphase-space distribution of dark matter around galaxyclusters [47, 63], for galaxy infall kinematics [71], as wellas in the spin-up of galactic halos [153] and proposed asuseful test of gravity with expected constraints of the or-der of | f R | . (10 − − − ). Particularly, the differencesbetween dynamically inferred masses and lensing massesof unshielded astronomical objects (see Sec. III D 4) andthe equivalence in the shielded case due to environmentalor self-shielding effects have been pointed out as promis-ing tests of gravity [27]. Zhao et al . [154] quantified theenvironmental dependence on the relation between thedynamical and lensing masses of dark matter halos, us-ing an indicator for the environment based on distancesto neighbouring halos and divided halos from N -bodysimulations into two different samples, which are eitherisolated or in high-density environments. They proposethat using their estimator and performing a similar di-vision in observed galaxy samples could yield a smoking8gun for modified gravity if a correlation between the en-vironment and the difference of dynamical and lensingmass is found.Dwarf galaxies in low-density environments are objectsof particular interest for the search of chameleon modifi-cations of gravity and to infer constraints on the model. Ifthe chameleon field amplitude is larger than the potentialwell of the dwarf galaxy and if its environmental densityis low enough that it does not shield it, the chameleonforce enhances the rotation curves of gas with respect tothe self-screened stars, yields a displacement between thetwo disk, and furthermore, warps the stellar disk and in-troduces an asymmetry in its rotation curve [41]. Thesesignatures can eventually be used to place constraints onthe chameleon modification of the order of the potentialwells of the dwarf galaxies, i.e., | f R | . − .Besides the gravitational modifications in the dwarfgalaxies that are expected to be strongest if the galax-ies reside in voids [27, 41], Martino and Sheth [155] alsopointed out an increase of the abundance of large voidsfor attractive extra forces, considering a Yukawa-typemodification of gravity. Li et al . [133] examined thisabundance in N -body simulations of f ( R ) gravity, con-firming the result of Ref. [155], also finding that halos areless screened near voids and that halos in voids are un-screened. Clampitt et al . [156] then analysed void statis-tics in chameleon gravity using excursion set theory topredict the increase in the number density of large voidsand furthermore, used this framework to explore the envi-ronmental dependence of void properties. Thereby, smallvoids in high-density environments were found to be emp-tier with faster expanding shells, motivating a cluster-ing analysis of small void tracers in redshift-space as adiscriminating test of gravity. In general, deviations invoid properties due to the modification of gravity arestronger than in halo properties, which could potentiallyyield more powerful tests of gravity with future data thanfrom constraining deviations in halo statistics.It is important to note that in order to exploit themodified gravity effects in these observables for tests ofgravity, in addition to running N -body simulations, ef-ficient modelling techniques such as based on excursionset theory or semi-analytic approaches as discussed inSec. III need to be further developed to allow for smoothvariations in the model parameter space and to prop-erly account for parameter degeneracies in the signatures.These modelling frameworks become even more substan-tial as the parameter space grows with the additionaldegrees of freedom introduced by the modifications ofgravity and emulator approaches become unfeasible.Ultimately, the techniques described in Sec. III maypotentially be generalised to describe the nonlinear large-scale structure formed in the full Horndeski theory de-scribed by Eq. (1). Important advances in this direc-tion have already been realised. The full linear cosmo-logical perturbations of the Horndeski theory have beenderived [157] and implemented in a Boltzmann lineartheory solver [158]. Note that the quasistatic pertur- bations take a relatively simple form and can easily bedetermined by the modifications of gravity in the back-ground. On nonlinear scales, in parallel to the descrip-tion of the cosmological small-scale structure formed inthe chameleon model, semi-analytic methods have alsobeen developed for other subclasses of the Horndeskiaction such as for the Galileon and symmetron mod-els [108, 159–161]. Barreira et al . [159] formulated aspherical collapse model for Galileon models, applyingtheir results in Ref. [160] to model the halo mass function,halo bias, and halo model matter power spectrum, andfinding good agreement with the corresponding statis-tics extracted from the N -body simulations performedin Ref. [162, 163]. Similarly, Schmidt et al . [164] stud-ied the spherical collapse model in the Dvali-Gabadadze-Porrati (DGP) model [165], which reduces to a Galileonmodel in the decoupling limit, and used it to determinehalo properties, comparing them to N -body simulationsconducted in Ref. [166]. Note that strong constraintson DGP and Galileon models have been inferred using avariety of linear cosmological observations in Refs. [167–170] and [160, 171, 172], respectively, limiting the impactof viable deviations from the concordance model on non-linear scales. Finally, the spherical collapse model, halomass function, halo bias, and nonlinear matter powerspectrum has also been analysed for symmetron mod-els [108, 161], a further subclass of the Horndeski action(see Sec. II).While the development of these modelling techniquesis still at an early stage, the success in reproducingthe important features of the modified gravity modelsin the nonlinear cosmological structure observed in N -body simulations of the models promises an applicationof these results not just to the Horndeski action. Itmay further provide hints for extending generalised linearcosmological perturbation formalisms such as the effec-tive field theory of cosmic acceleration [157, 173] or thePPF framework [174], which both embed the linear per-turbations of the Horndeski theory, to nonlinear scales.Similarly, the success of the description of the nonlinearstructure in Galileon and DGP gravity based on thesesemi-analytic models may also be shared by massive grav-ity [175], providing an interesting area of study for futurework. VI. CONCLUSION
The presence of a chameleon field in our Universe maymodify the gravitational interactions on cluster scaleswhile recovering general relativity locally. Whereas N -body simulations provide an essential tool for studyingthe effects of the scalar-tensor modification of gravity, es-pecially in regions where the modification is getting sup-pressed due to the chameleon mechanism, semi-analyticmodelling techniques become a necessity for the compar-ison of the theoretical signatures to observational data.These tests require an efficient interpolation and extrapo-9lation of the simulated results between different choices ofcosmological parameters and model specifications of thegravitational theories of interest. Emulator approacheswill not be feasible for this task as the parameter spacegrows with the extra degrees of freedom introduced bythe different modified gravity models and the computa-tional methods need to be generalised accordingly.This article summarises a range of observable signa-tures in the nonlinear cosmological structure that arecharacteristic for the chameleon modification. It reviewsand compares different techniques to model these observ-ables based on N -body simulations, different phenomeno-logical formalisms, fitting functions to simulations, ana-lytical and numerical approximations, the spherical col-lapse model, excursion set theory, the halo model, andperturbation theory. Hereby, a particular focus is givento the well studied Hu-Sawicki and designer models of f ( R ) gravity, for which a summary of the current state ofobservational constraints is provided and supplementedwith an outlook on prospective constraints and novelmethods for testing the chameleon modifications.The semi-analytic methods discussed here are still atan early stage of construction and need to be developed further. Their success in reproducing the important fea-tures of the nonlinear cosmological structure observedin N -body simulations of chameleon f ( R ) gravity, alongwith similar achievements for describing the cosmologi-cal structure in Galileon and symmetron models, antici-pates that a generalisation of these modelling techniquesand an application thereof to the full Horndeski theory ofscalar-tensor gravity may be feasible. Hence, at the near-ing of the 100th anniversary since the formulation of thefoundations of general relativity [176], the study of thecosmological structure formed in the simplest extensionof the Theory of Gravity, i.e., for scalar-tensor models,and the observational constraints that can be inferred onthese extensions promise to remain a very interesting andactive field of research. a. Acknowledgements The author thanks PhilippeBrax, Baojiu Li, Patrik Valageas, and Gong-Bo Zhao forsharing numerical results which have been used to pro-duce Figs. 3, 4, and 5. This work has been supportedby the STFC Consolidated Grant for Astronomy andAstrophysics of the University of Edinburgh. Numer-ical computations have been performed with Wolfram
M athematica
9. Please contact the author for accessto research materials. [1] C. M. Will,
The confrontation between generalrelativity and experiment , Living Rev. Rel. (2005) 3,[ gr-qc/0510072 ].[2] J. Khoury and A. Weltman, Chameleon fields:Awaiting surprises for tests of gravity in space , Phys.Rev.Lett. (2004) 171104, [ astro-ph/0309300 ].[3] J. Khoury and A. Weltman, Chameleon cosmology , Phys. Rev.
D69 (2004) 044026, [ astro-ph/0309411 ].[4] P. Brax, C. van de Bruck, A.-C. Davis, J. Khoury, andA. Weltman,
Detecting dark energy in orbit - TheCosmological chameleon , Phys.Rev.
D70 (2004)123518, [ astro-ph/0408415 ].[5] J. Cembranos,
The Newtonian limit at intermediateenergies , Phys.Rev.
D73 (2006) 064029,[ gr-qc/0507039 ].[6] D. F. Mota and D. J. Shaw,
Evading EquivalencePrinciple Violations, Cosmological and otherExperimental Constraints in Scalar Field Theories witha Strong Coupling to Matter , Phys.Rev.
D75 (2007)063501, [ hep-ph/0608078 ].[7] J. Wang, L. Hui, and J. Khoury,
No-Go Theorems forGeneralized Chameleon Field Theories , Phys.Rev.Lett. (2012) 241301, [ arXiv:1208.4612 ].[8] W. Hu and I. Sawicki,
Models of f(R) CosmicAcceleration that Evade Solar-System Tests , Phys. Rev.
D76 (2007) 064004, [ arXiv:0705.1158 ].[9] Y.-S. Song, W. Hu, and I. Sawicki,
The large scalestructure of f(R) gravity , Phys. Rev.
D75 (2007)044004, [ astro-ph/0610532 ].[10] H. A. Buchdahl,
Non-linear Lagrangians andcosmological theory , Mon. Not. Roy. Astron. Soc. (1970) 1. [11] A. A. Starobinsky,
Spectrum of relict gravitationalradiation and the early state of the universe , JETPLett. (1979) 682–685.[12] A. A. Starobinsky, A new type of isotropic cosmologicalmodels without singularity , Phys. Lett.
B91 (1980)99–102.[13] S. Capozziello, S. Carloni, and A. Troisi,
Quintessencewithout scalar fields , Recent Res. Dev. Astron.Astrophys. (2003) 625, [ astro-ph/0303041 ].[14] S. M. Carroll, V. Duvvuri, M. Trodden, and M. S.Turner, Is cosmic speed-up due to new gravitationalphysics? , Phys. Rev.
D70 (2004) 043528,[ astro-ph/0306438 ].[15] S. Nojiri and S. D. Odintsov,
Modified gravity withnegative and positive powers of the curvature:Unification of the inflation and of the cosmicacceleration , Phys. Rev.
D68 (2003) 123512,[ hep-th/0307288 ].[16] T. Chiba, , Phys.Lett.
B575 (2003) 1–3, [ astro-ph/0307338 ].[17] P. Zhang,
Testing f ( R ) gravity against the large scalestructure of the universe. , Phys.Rev.
D73 (2006)123504, [ astro-ph/0511218 ].[18] L. Amendola, D. Polarski, and S. Tsujikawa,
Are f(R)dark energy models cosmologically viable ? , Phys.Rev.Lett. (2007) 131302, [ astro-ph/0603703 ].[19] L. Amendola, D. Polarski, and S. Tsujikawa, Power-laws f(R) theories are cosmologicallyunacceptable , Int.J.Mod.Phys.
D16 (2007) 1555–1561,[ astro-ph/0605384 ].[20] L. Amendola, R. Gannouji, D. Polarski, andS. Tsujikawa,
Conditions for the cosmological viabilityof f(R) dark energy models , Phys.Rev.
D75 (2007) gr-qc/0612180 ].[21] B. Li and J. D. Barrow, The Cosmology of f(R) gravityin metric variational approach , Phys.Rev.
D75 (2007)084010, [ gr-qc/0701111 ].[22] P. Zhang, M. Liguori, R. Bean, and S. Dodelson,
Probing Gravity at Cosmological Scales byMeasurements which Test the Relationship betweenGravitational Lensing and Matter Overdensity , Phys.Rev. Lett. (2007) 141302, [ arXiv:0704.1932 ].[23] Y.-S. Song, H. Peiris, and W. Hu, CosmologicalConstraints on f(R) Acceleration Models , Phys. Rev.
D76 (2007) 063517, [ arXiv:0706.2399 ].[24] B. Jain and P. Zhang,
Observational Tests of ModifiedGravity , Phys.Rev.
D78 (2008) 063503,[ arXiv:0709.2375 ].[25] P. Brax, C. van de Bruck, A.-C. Davis, and D. J.Shaw, f(R) Gravity and Chameleon Theories , Phys.Rev.
D78 (2008) 104021, [ arXiv:0806.3415 ].[26] G.-B. Zhao, L. Pogosian, A. Silvestri, andJ. Zylberberg,
Searching for modified growth patternswith tomographic surveys , Phys.Rev.
D79 (2009)083513, [ arXiv:0809.3791 ].[27] L. Hui, A. Nicolis, and C. Stubbs,
EquivalencePrinciple Implications of Modified Gravity Models , Phys.Rev.
D80 (2009) 104002, [ arXiv:0905.2966 ].[28] T. L. Smith,
Testing gravity on kiloparsec scales withstrong gravitational lenses , arXiv:0907.4829 .[29] F. Schmidt, A. Vikhlinin, and W. Hu, ClusterConstraints on f(R) Gravity , Phys. Rev.
D80 (2009)083505, [ arXiv:0908.2457 ].[30] T. Giannantonio, M. Martinelli, A. Silvestri, andA. Melchiorri,
New constraints on parametrisedmodified gravity from correlations of the CMB withlarge scale structure , JCAP (2010) 030,[ arXiv:0909.2045 ].[31] R. Reyes et al.,
Confirmation of general relativity onlarge scales from weak lensing and galaxy velocities , Nature (2010.) 256–258, [ arXiv:1003.2185 ].[32] L. Lombriser, A. Slosar, U. Seljak, and W. Hu,
Constraints on f(R) gravity from probing thelarge-scale structure , Phys. Rev.
D85 (2012) 124038,[ arXiv:1003.3009 ].[33] K. Yamamoto, G. Nakamura, G. Hutsi, T. Narikawa,and T. Sato,
Constraint on the cosmological f(R)model from the multipole power spectrum of the SDSSluminous red galaxy sample and prospects for a futureredshift survey , Phys.Rev.
D81 (2010) 103517,[ arXiv:1004.3231 ].[34] B. Jain and J. Khoury,
Cosmological Tests of Gravity , Annals Phys. (2010) 1479–1516,[ arXiv:1004.3294 ].[35] H. Motohashi, A. A. Starobinsky, and J. Yokoyama,
Matter power spectrum in f(R) gravity with massiveneutrinos , Prog.Theor.Phys. (2010) 541–546,[ arXiv:1005.1171 ].[36] X. Gao,
Testing gravity with non-Gaussianity , Phys.Lett.
B702 (2011) 197–200, [ arXiv:1008.2123 ].[37] X. Wang, X. Chen, and C. Park,
Topology of largescale structure as test of modified gravity , Astrophys.J. (2012) 48, [ arXiv:1010.3035 ].[38] S. Ferraro, F. Schmidt, and W. Hu,
Cluster Abundancein f(R) Gravity Models , Phys.Rev.
D83 (2011) 063503,[ arXiv:1011.0992 ]. [39] S. A. Thomas, S. A. Appleby, and J. Weller,
ModifiedGravity: the CMB, Weak Lensing and GeneralParameterisations , JCAP (2011) 036,[ arXiv:1101.0295 ].[40] A.-C. Davis, E. A. Lim, J. Sakstein, and D. Shaw,
Modified Gravity Makes Galaxies Brighter , Phys.Rev.
D85 (2012) 123006, [ arXiv:1102.5278 ].[41] B. Jain and J. VanderPlas,
Tests of Modified Gravitywith Dwarf Galaxies , JCAP (2011) 032,[ arXiv:1106.0065 ].[42] A. Hojjati, L. Pogosian, and G.-B. Zhao,
Testinggravity with CAMB and CosmoMC , JCAP (2011) 005, [ arXiv:1106.4543 ].[43] G.-B. Zhao, H. Li, E. V. Linder, K. Koyama, D. J.Bacon, et al.,
Testing Einstein Gravity with CosmicGrowth and Expansion , Phys.Rev.
D85 (2012) 123546,[ arXiv:1109.1846 ].[44] H. Gil-Marin, F. Schmidt, W. Hu, R. Jimenez, andL. Verde,
The Bispectrum of f(R) Cosmologies , JCAP (2011) 019, [ arXiv:1109.2115 ].[45] R. Wojtak, S. H. Hansen, and J. Hjorth,
Gravitationalredshift of galaxies in clusters as predicted by generalrelativity , Nature (2011) 567–569,[ arXiv:1109.6571 ].[46] L. Lombriser, F. Schmidt, T. Baldauf,R. Mandelbaum, U. Seljak, et al.,
Cluster DensityProfiles as a Test of Modified Gravity , Phys. Rev.
D85 (2012) 102001, [ arXiv:1111.2020 ].[47] T. Y. Lam, T. Nishimichi, F. Schmidt, and M. Takada,
Testing Gravity with the Stacked Phase Space aroundGalaxy Clusters , Phys.Rev.Lett. (2012) 051301,[ arXiv:1202.4501 ].[48] A. Terukina and K. Yamamoto,
Gas density profile indark matter halo in chameleon cosmology , Phys.Rev.
D86 (2012) 103503, [ arXiv:1203.6163 ].[49] E. Di Valentino, A. Melchiorri, V. Salvatelli, andA. Silvestri,
Parametrised modified gravity and theCMB Bispectrum , Phys.Rev.
D86 (2012) 063517,[ arXiv:1204.5352 ].[50] B. Jain, V. Vikram, and J. Sakstein,
AstrophysicalTests of Modified Gravity: Constraints from DistanceIndicators in the Nearby Universe , Astrophys.J. (2013) 39, [ arXiv:1204.6044 ].[51] L. Samushia, B. A. Reid, M. White, W. J. Percival,A. J. Cuesta, et al.,
The Clustering of Galaxies in theSDSS-III DR9 Baryon Oscillation SpectroscopicSurvey: Testing Deviations from Λ and GeneralRelativity using anisotropic clustering of galaxies , Mon. Not. Roy. Astron. Soc. (2013) 1514–1528,[ arXiv:1206.5309 ].[52] J.-h. He,
Testing f ( R ) dark energy model with thelarge scale structure , Phys.Rev.
D86 (2012) 103505,[ arXiv:1207.4898 ].[53] B. Hu, M. Liguori, N. Bartolo, and S. Matarrese,
Future CMB ISW-Lensing bispectrum constraints onmodified gravity in the Parameterized Post-Friedmannformalism , Phys.Rev.
D88 (2013) 024012,[ arXiv:1211.5032 ].[54] H. Okada, T. Totani, and S. Tsujikawa,
Constraints onf(R) theory and Galileons from the latest data ofgalaxy redshift surveys , Phys.Rev.
D87 (2013) 103002,[ arXiv:1208.4681 ].[55] A. Hall, C. Bonvin, and A. Challinor,
Testing GeneralRelativity with 21 cm intensity mapping , Phys.Rev. D87 (2013) 064026, [ arXiv:1212.0728 ].[56] F. Simpson, C. Heymans, D. Parkinson, C. Blake,M. Kilbinger, et al.,
CFHTLenS: testing the laws ofgravity with tomographic weak lensing andredshift-space distortions , Mon. Not. Roy. Astron. Soc. (Mar., 2013) 2249–2263, [ arXiv:1212.3339 ].[57] F. Albareti, J. Cembranos, A. de la Cruz-Dombriz,and A. Dobado,
On the non-attractive character ofgravity in f(R) theories , JCAP (2013) 009,[ arXiv:1212.4781 ].[58] L. Lombriser, J. Yoo, and K. Koyama,
Relativisticeffects in galaxy clustering in a parametrizedpost-Friedmann universe , Phys.Rev.
D87 (2013)104019, [ arXiv:1301.3132 ].[59] P. Brax, A.-C. Davis, and J. Sakstein,
PulsarConstraints on Screened Modified Gravity , arXiv:1301.5587 .[60] A. Marchini, A. Melchiorri, V. Salvatelli, andL. Pagano, Constraints on Modified Gravity from ACTand SPT , Phys.Rev.
D87 (2013) 083527,[ arXiv:1302.2593 ].[61] V. Vikram, A. Cabr´e, B. Jain, and J. VanderPlas,
Astrophysical Tests of Modified Gravity: theMorphology and Kinematics of Dwarf Galaxies , JCAP (2013) 020, [ arXiv:1303.0295 ].[62] A. Abebe, A. de la Cruz-Dombriz, and P. K. S.Dunsby,
Large Scale Structure Constraints for a Classof f(R) Theories of Gravity , Phys.Rev.
D88 (2013)044050, [ arXiv:1304.3462 ].[63] T. Y. Lam, F. Schmidt, T. Nishimichi, and M. Takada,
Modeling the Phase-Space Distribution around MassiveHalos , Phys.Rev.
D88 (2013) 023012,[ arXiv:1305.5548 ].[64] W. A. Hellwing, B. Li, C. S. Frenk, and S. Cole,
Hierarchical clustering in chameleon f(R) gravity , Mon. Not. Roy. Astron. Soc. (Nov., 2013)2806–2821, [ arXiv:1305.7486 ].[65] A. Upadhye and J. H. Steffen,
Monopole radiation inmodified gravity , arXiv:1306.6113 .[66] A. Marchini and V. Salvatelli, Updated constraintsfrom the PLANCK experiment on modified gravity , Phys.Rev.
D88 (2013) 027502, [ arXiv:1307.2002 ].[67] J.-h. He,
Weighing Neutrinos in f ( R ) gravity , Phys.Rev. D 88, (2013) 103523, [ arXiv:1307.4876 ].[68] B. Hu, M. Liguori, N. Bartolo, and S. Matarrese,
Parametrized modified gravity constraints after Planck , Phys.Rev.
D88 (2013) 123514, [ arXiv:1307.5276 ].[69] N. Mirzatuny, S. Khosravi, S. Baghram, andH. Moshafi,
Simultaneous effect of modified gravity andprimordial non-Gaussianity in large scale structureobservations , JCAP (2014) 019,[ arXiv:1308.2874 ].[70] J. Sakstein,
Stellar Oscillations in Modified Gravity , Phys.Rev.
D88 (2013) 124013, [ arXiv:1309.0495 ].[71] Y. Zu, D. Weinberg, E. Jennings, B. Li, andM. Wyman,
Galaxy Infall Kinematics as a Test ofModified Gravity , arXiv:1310.6768 .[72] Y.-C. Cai, B. Li, S. Cole, C. S. Frenk, andM. Neyrinck, The integrated Sachs-Wolfe effect in f(R)gravity , Mon. Not. Roy. Astron. Soc. (Feb., 2014)[ arXiv:1310.6986 ].[73] M. Baldi, F. Villaescusa-Navarro, M. Viel,E. Puchwein, V. Springel, et al.,
Cosmic DegeneraciesI: Joint N-body Simulations of Modified Gravity and Massive Neutrinos , arXiv:1311.2588 .[74] Heidelberg Institute for Theoretical Studies
Collaboration,
Scaling relations and mass bias inhydrodynamical f(R) gravity simulations of galaxyclusters , arXiv:1311.5560 .[75] L. Lombriser, K. Koyama, and B. Li, Halo modellingin chameleon theories , JCAP (2014) 021,[ arXiv:1312.1292 ].[76] P. Brax and A. Upadhye,
Chameleon Fragmentation , JCAP (2014) 018, [ arXiv:1312.2747 ].[77] A. Terukina, L. Lombriser, K. Yamamoto, D. Bacon,K. Koyama, et al.,
Testing chameleon gravity with theComa cluster , arXiv:1312.5083 .[78] W. A. Hellwing, A. Barreira, C. S. Frenk, B. Li, andS. Cole, A clear and measurable signature of modifiedgravity in the galaxy velocity field , arXiv:1401.0706 .[79] J. Dossett, B. Hu, and D. Parkinson, Constrainingmodels of f ( R ) gravity with Planck and WiggleZ powerspectrum data , arXiv:1401.3980 .[80] A.-C. Davis, R. Gregory, R. Jha, and J. Muir, Astrophysical black holes in screened modified gravity , arXiv:1402.4737 .[81] D. Munshi, B. Hu, A. Renzi, A. Heavens, and P. Coles, Probing Modified Gravity Theories with ISW and CMBLensing , arXiv:1403.0852 .[82] A. L. Erickcek, N. Barnaby, C. Burrage, and Z. Huang, Catastrophic Consequences of Kicking the Chameleon , Phys. Rev. Lett. (2013)[ arXiv:1304.0009 ].[83] A. L. Erickcek, N. Barnaby, C. Burrage, and Z. Huang,
Chameleons in the Early Universe: Kicks, Rebounds,and Particle Production , arXiv:1310.5149 .[84] A. Upadhye, W. Hu, and J. Khoury, Quantum Stabilityof Chameleon Field Theories , Phys.Rev.Lett. (2012) 041301, [ arXiv:1204.3906 ].[85] H. Oyaizu,
Non-linear evolution of f(R) cosmologies I:methodology , Phys. Rev.
D78 (2008) 123523,[ arXiv:0807.2449 ].[86] H. Oyaizu, M. Lima, and W. Hu,
Non-linear evolutionof f(R) cosmologies II: power spectrum , Phys. Rev.
D78 (2008) 123524, [ arXiv:0807.2462 ].[87] B. Li and H. Zhao,
Structure Formation by Fifth ForceI: N-Body vs. Linear Simulations , Phys.Rev.
D80 (2009) 044027, [ arXiv:0906.3880 ].[88] B. Li, D. F. Mota, and J. D. Barrow,
N-bodySimulations for Extended Quintessence Models , Astrophys.J. (2011) 109, [ arXiv:1009.1400 ].[89] G.-B. Zhao, B. Li, and K. Koyama,
N-bodySimulations for f(R) Gravity using a Self-adaptiveParticle-Mesh Code , Phys. Rev.
D83 (2011) 044007,[ arXiv:1011.1257 ].[90] B. Li, G.-B. Zhao, R. Teyssier, and K. Koyama,
ECOSMOG: An Efficient Code for Simulating ModifiedGravity , JCAP (2012) 051, [ arXiv:1110.1379 ].[91] E. Puchwein, M. Baldi, and V. Springel,
Modified-Gravity-GADGET: a new code forcosmological hydrodynamical simulations of modifiedgravity models , Mon. Not. Roy. Astron. Soc. (Nov., 2013) 348–360, [ arXiv:1305.2418 ].[92] C. Llinares, D. F. Mota, and H. A. Winther,
ISIS: anew N-body cosmological code with scalar fields basedon RAMSES , arXiv:1307.6748 .[93] W. Hu and I. Sawicki, A ParameterizedPost-Friedmann Framework for Modified Gravity , Phys. Rev.
D76 (2007) 104043, [ arXiv:0708.1190 ].[94] G.-B. Zhao et al.,
Probing modifications of GeneralRelativity using current cosmological observations , Phys. Rev.
D81 (2010) 103510, [ arXiv:1003.0001 ].[95] Y. Li and W. Hu,
Chameleon Halo Modeling in f(R)Gravity , Phys.Rev.
D84 (2011) 084033,[ arXiv:1107.5120 ].[96] G.-B. Zhao,
Modeling the nonlinear clustering inmodified gravity models I: A fitting formula for matterpower spectrum of f(R) gravity , arXiv:1312.1291 .[97] F. Schmidt, Dynamical Masses in Modified Gravity , Phys. Rev.
D81 (2010) 103002, [ arXiv:1003.0409 ].[98] R. Pourhasan, N. Afshordi, R. Mann, and A. Davis,
Chameleon Gravity, Electrostatics, and Kinematics inthe Outer Galaxy , JCAP (2011) 005,[ arXiv:1109.0538 ].[99] L. Lombriser, K. Koyama, G.-B. Zhao, and B. Li,
Chameleon f(R) gravity in the virialized cluster , Phys.Rev.
D85 (2012) 124054, [ arXiv:1203.5125 ].[100] F. Schmidt, M. V. Lima, H. Oyaizu, and W. Hu,
Non-linear Evolution of f(R) Cosmologies III: HaloStatistics , Phys. Rev.
D79 (2009) 083518,[ arXiv:0812.0545 ].[101] A. Borisov, B. Jain, and P. Zhang,
Spherical Collapsein f(R) Gravity , Phys.Rev.
D85 (2012) 063518,[ arXiv:1102.4839 ].[102] B. Li and G. Efstathiou,
An Extended Excursion SetApproach to Structure Formation in ChameleonModels , Mon. Not. Roy. Astron. Soc. (2012) 1431,[ arXiv:1110.6440 ].[103] L. Lombriser, B. Li, K. Koyama, and G.-B. Zhao,
Modeling halo mass functions in chameleon f(R)gravity , Phys. Rev. D 87, (2013)[ arXiv:1304.6395 ].[104] B. Li and T. Y. Lam,
Excursion set theory for modifiedgravity: Eulerian versus Lagrangian environments , Mon. Not. Roy. Astron. Soc. (Sept., 2012)730–739, [ arXiv:1205.0058 ].[105] T. Y. Lam and B. Li,
Excursion set theory for modifiedgravity: correlated steps, mass functions and halo bias , Mon. Not. Roy. Astron. Soc. (Nov., 2012)3260–3270, [ arXiv:1205.0059 ].[106] M. Kopp, S. A. Appleby, I. Achitouv, and J. Weller,
Spherical collapse and halo mass function in f(R)theories , Phys.Rev.
D88 (2013) 084015,[ arXiv:1306.3233 ].[107] K. Koyama, A. Taruya, and T. Hiramatsu,
Non-linearEvolution of Matter Power Spectrum in ModifiedTheory of Gravity , Phys.Rev.
D79 (2009) 123512,[ arXiv:0902.0618 ].[108] P. Brax and P. Valageas,
Impact on the powerspectrum of screening in modified gravity scenarios , Phys.Rev. (July, 2013) 023527, [ arXiv:1305.5647 ].[109] G. W. Horndeski, Second-order scalar-tensor fieldequations in a four-dimensional space , Int.J.Theor.Phys. (1974) 363–384.[110] C. Deffayet, X. Gao, D. Steer, and G. Zahariade, Fromk-essence to generalised Galileons , Phys.Rev.
D84 (2011) 064039, [ arXiv:1103.3260 ].[111] A. Nicolis, R. Rattazzi, and E. Trincherini,
TheGalileon as a local modification of gravity , Phys.Rev.
D79 (2009) 064036, [ arXiv:0811.2197 ].[112] K. Hinterbichler and J. Khoury,
Symmetron Fields:Screening Long-Range Forces Through Local Symmetry Restoration , Phys. Rev. Lett. (2010) 231301,[ arXiv:1001.4525 ].[113] A. Upadhye,
Dark energy fifth forces in torsionpendulum experiments , Phys.Rev.
D86 (2012) 102003,[ arXiv:1209.0211 ].[114] L. Pogosian and A. Silvestri,
The pattern of growth inviable f(R) cosmologies , Phys. Rev.
D77 (2008)023503, [ arXiv:0709.0296 ].[115] S. Nojiri and S. D. Odintsov,
Modified f(R) gravityconsistent with realistic cosmology: From matterdominated epoch to dark energy universe , Phys.Rev.
D74 (2006) 086005, [ hep-th/0608008 ].[116] S. Nojiri and S. D. Odintsov,
Modified gravity and itsreconstruction from the universe expansion history , J.Phys.Conf.Ser. (2007) 012005, [ hep-th/0611071 ].[117] S. Nojiri and S. D. Odintsov, Unified cosmic history inmodified gravity: from F(R) theory to Lorentznon-invariant models , Phys.Rept. (2011) 59–144,[ arXiv:1011.0544 ].[118] G. Esposito-Farese and D. Polarski,
Scalar tensorgravity in an accelerating universe , Phys.Rev.
D63 (2001) 063504, [ gr-qc/0009034 ].[119] S. Tsujikawa, K. Uddin, S. Mizuno, R. Tavakol, andJ. Yokoyama,
Constraints on scalar-tensor models ofdark energy from observational and local gravity tests , Phys.Rev.
D77 (2008) 103009, [ arXiv:0803.1106 ].[120] A. Hojjati, L. Pogosian, A. Silvestri, and S. Talbot,
Practical solutions for perturbed f(R) gravity , Phys.Rev.
D86 (2012) 123503, [ arXiv:1210.6880 ].[121] N. A. Lima and A. R. Liddle,
Linear perturbations inviable f(R) theories , Phys. Rev. D 88, (2013)[ arXiv:1307.1613 ].[122] J. Noller, F. von Braun-Bates, and P. G. Ferreira,
Relativistic scalar fields and the quasi-staticapproximation in theories of modified gravity , Phys.Rev.
D89 (2014) 023521, [ arXiv:1310.3266 ].[123] L. Lombriser,
Consistency check of Λ CDMphenomenology , Phys. Rev.
D83 (2011) 063519,[ arXiv:1101.0594 ].[124] A. Knebe, A. Green, and J. Binney,
Mlapm - a c codefor cosmological simulations , Mon. Not. Roy. Astron.Soc. (2001) 845, [ astro-ph/0103503 ].[125] B. Li, W. A. Hellwing, K. Koyama, G.-B. Zhao,E. Jennings, et al.,
The nonlinear matter and velocitypower spectra in f(R) gravity , Mon. Not. Roy. Astron.Soc. (2013) 743–755, [ arXiv:1206.4317 ].[126] R. Teyssier,
Cosmological hydrodynamics with adaptivemesh refinement: a new high resolution code calledramses , Astron.Astrophys. (2002) 337–364,[ astro-ph/0111367 ].[127] V. Springel,
The Cosmological simulation codeGADGET-2 , Mon. Not. Roy. Astron. Soc. (2005)1105–1134, [ astro-ph/0505010 ].[128] C. Llinares and D. Mota,
Releasing scalar fields:cosmological simulations of scalar-tensor theories forgravity beyond the static approximation , Phys.Rev.Lett. (2013), no. 16 161101, [ arXiv:1302.1774 ].[129] E. Jennings, C. M. Baugh, B. Li, G.-B. Zhao, andK. Koyama,
Redshift space distortions in f(R) gravity , Mon. Not. Roy. Astron. Soc. (2012) 2128–2143,[ arXiv:1205.2698 ].[130] M. B. Gronke, C. Llinares, and D. F. Mota,
Gravitational redshift profiles in the f(R) andsymmetron models , Astron. Astrophys. (Feb., arXiv:1307.6994 ].[131] F. Fontanot, E. Puchwein, V. Springel, andD. Bianchi, Semi-analytic galaxy formation inf(R)-gravity cosmologies , Mon. Not. Roy. Astron. Soc. (Dec., 2013) 2672–2679, [ arXiv:1307.5065 ].[132] P. Brax, R. Rosenfeld, and D. Steer,
SphericalCollapse in Chameleon Models , JCAP (2010)033, [ arXiv:1005.2051 ].[133] B. Li, G.-B. Zhao, and K. Koyama,
Haloes and voidsin f(R) gravity , Mon. Not. Roy. Astron. Soc. (Apr., 2012) 3481–3487, [ arXiv:1111.2602 ].[134] T. Y. Lam and R. K. Sheth,
Perturbation theory andexcursion set estimates of the probability distributionfunction of dark matter, and a method forreconstructing the initial distribution function , Mon.Not. Roy. Astron. Soc. (May, 2008) 407–415,[ arXiv:0711.5029 ].[135] W. H. Press and P. Schechter,
Formation of Galaxiesand Clusters of Galaxies by Self-Similar GravitationalCondensation , Astrophys. J. (Feb., 1974) 425–438.[136] R. K. Sheth and G. Tormen,
Large scale bias and thepeak background split , Mon. Not. Roy. Astron. Soc. (1999) 119, [ astro-ph/9901122 ].[137] R. K. Sheth, H. Mo, and G. Tormen,
Ellipsoidalcollapse and an improved model for the number andspatial distribution of dark matter haloes , Mon. Not.Roy. Astron. Soc. (2001) 1, [ astro-ph/9907024 ].[138] R. K. Sheth and G. Tormen,
An Excursion set modelof hierarchical clustering : Ellipsoidal collapse and themoving barrier , Mon. Not. Roy. Astron. Soc. (2002) 61, [ astro-ph/0105113 ].[139] L. Hui and K. P. Parfrey,
The Evolution of Bias:Generalized , Phys.Rev.
D77 (2008) 043527,[ arXiv:0712.1162 ].[140] M. C. Martino, H. F. Stabenau, and R. K. Sheth,
Spherical Collapse and Modified Gravity , Phys. Rev.
D79 (2009) 084013, [ arXiv:0812.0200 ].[141] K. Parfrey, L. Hui, and R. K. Sheth,
Scale-dependenthalo bias from scale-dependent growth , Phys.Rev.
D83 (2011) 063511, [ arXiv:1012.1335 ].[142] J. F. Navarro, C. S. Frenk, and S. D. M. White,
TheStructure of Cold Dark Matter Halos , Astrophys. J. (1996) 563–575, [ astro-ph/9508025 ].[143] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville,A. V. Kravtsov, et al.,
Profiles of dark haloes.Evolution, scatter, and environment , Mon. Not. Roy.Astron. Soc. (2001) 559–575, [ astro-ph/9908159 ].[144] J. Peacock and R. Smith,
Halo occupation numbersand galaxy bias , Mon. Not. Roy. Astron. Soc. (2000) 1144, [ astro-ph/0005010 ].[145] U. Seljak,
Analytic model for galaxy and dark matterclustering , Mon. Not. Roy. Astron. Soc. (2000)203, [ astro-ph/0001493 ].[146] A. Cooray and R. K. Sheth,
Halo models of large scalestructure , Phys. Rept. (2002) 1–129,[ astro-ph/0206508 ].[147] D. J. Eisenstein and W. Hu,
Baryonic features in thematter transfer function , Astrophys.J. (1998) 605,[ astro-ph/9709112 ].[148] D. J. Eisenstein and W. Hu,
Power spectra for colddark matter and its variants , Astrophys.J. (1997)5, [ astro-ph/9710252 ].[149]
Virgo Consortium
Collaboration, R. Smith et al.,
Stable clustering, the halo model and nonlinear cosmological power spectra , Mon. Not. Roy. Astron.Soc. (2003) 1311, [ astro-ph/0207664 ].[150] P. Valageas, T. Nishimichi, and A. Taruya,
Matterpower spectrum from a Lagrangian-space regularizationof perturbation theory , Phys.Rev.
D87 (2013) 083522,[ arXiv:1302.4533 ].[151] C. Llinares and D. F. Mota,
Shape of Clusters ofGalaxies as a Probe of Screening Mechanisms inModified Gravity , Phys.Rev.Lett. (2013), no. 15151104, [ arXiv:1205.5775 ].[152] A. S. Bolton, S. Rappaport, and S. Burles,
Constrainton the Post-Newtonian Parameter gamma on GalacticSize Scales , Phys.Rev.
D74 (2006) 061501,[ astro-ph/0607657 ].[153] J. Lee, G.-B. Zhao, B. Li, and K. Koyama,
ModifiedGravity Spins Up Galactic Halos , Astrophys.J. (2013) 28, [ arXiv:1204.6608 ].[154] G.-B. Zhao, B. Li, and K. Koyama,
Testing GeneralRelativity using the Environmental Dependence ofDark Matter Halos , Phys.Rev.Lett. (2011) 071303,[ arXiv:1105.0922 ].[155] M. C. Martino and R. K. Sheth,
Density profiles andvoids in modified gravity models , arXiv:0911.1829 .[156] J. Clampitt, Y.-C. Cai, and B. Li, Voids in modifiedgravity: excursion set predictions , Mon. Not. Roy.Astron. Soc. (May, 2013) 749–766,[ arXiv:1212.2216 ].[157] A. De Felice, T. Kobayashi, and S. Tsujikawa,
Effective gravitational couplings for cosmologicalperturbations in the most general scalar-tensor theorieswith second-order field equations , Phys.Lett.
B706 (2011) 123–133, [ arXiv:1108.4242 ].[158] B. Hu, M. Raveri, N. Frusciante, and A. Silvestri,
Effective Field Theory of Cosmic Acceleration: animplementation in CAMB , arXiv:1312.5742 .[159] A. Barreira, B. Li, C. M. Baugh, and S. Pascoli, Spherical collapse in Galileon gravity: fifth forcesolutions, halo mass function and halo bias , JCAP (2013) 056, [ arXiv:1308.3699 ].[160] A. Barreira, B. Li, W. A. Hellwing, L. Lombriser,C. M. Baugh, et al.,
Halo model and halo properties inGalileon gravity cosmologies , arXiv:1401.1497 .[161] L. Taddei, R. Catena, and M. Pietroni, Sphericalcollapse and halo mass function in the symmetronmodel , arXiv:1310.6175 .[162] A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh, andS. Pascoli, Nonlinear structure formation in the cubicGalileon gravity model , JCAP (Oct., 2013) 27,[ arXiv:1306.3219 ].[163] B. Li, A. Barreira, C. M. Baugh, W. A. Hellwing,K. Koyama, et al., Simulating the quartic Galileongravity model on adaptively refined meshes , JCAP (2013) 012, [ arXiv:1308.3491 ].[164] F. Schmidt, W. Hu, and M. Lima,
Spherical Collapseand the Halo Model in Braneworld Gravity , Phys.Rev.
D81 (2010) 063005, [ arXiv:0911.5178 ].[165] G. Dvali, G. Gabadadze, and M. Porrati, , Phys.Lett.
B485 (2000) 208–214, [ hep-th/0005016 ].[166] F. Schmidt,
Cosmological Simulations ofNormal-Branch Braneworld Gravity , Phys.Rev.
D80 (2009) 123003, [ arXiv:0910.0235 ].[167] W. Fang et al.,
Challenges to the DGP Model fromHorizon-Scale Growth and Geometry , Phys. Rev.
D78 (2008) 103509, [ arXiv:0808.2208 ].[168] L. Lombriser, W. Hu, W. Fang, and U. Seljak, Cosmological Constraints on DGP Braneworld Gravitywith Brane Tension , Phys. Rev.
D80 (2009) 063536,[ arXiv:0905.1112 ].[169] A. Raccanelli, D. Bertacca, D. Pietrobon, F. Schmidt,L. Samushia, et al.,
Testing gravity using large-scaleredshift-space distortions , Mon. Not. Roy. Astron. Soc. (Nov., 2013) 89–100, [ arXiv:1207.0500 ].[170] L. Xu,
Confronting DGP braneworld gravity withcosmic observations after Planck data , JCAP (2014) 048, [ arXiv:1312.4679 ].[171] S. A. Appleby and E. V. Linder,
Trial of Galileongravity by cosmological expansion and growthobservations , JCAP (2012) 026,[ arXiv:1204.4314 ].[172] A. Barreira, B. Li, A. Sanchez, C. M. Baugh, andS. Pascoli,
The parameter space in Galileon gravity models , Phys.Rev.
D87 (2013) 103511,[ arXiv:1302.6241 ].[173] J. K. Bloomfield, E. E. Flanagan, M. Park, andS. Watson,
Dark energy or modified gravity? Aneffective field theory approach , JCAP (2013) 010,[ arXiv:1211.7054 ].[174] T. Baker, P. G. Ferreira, and C. Skordis,
TheParameterized Post-Friedmann Framework forTheories of Modified Gravity: Concepts, Formalismand Examples , Phys.Rev.
D87 (2013) 024015,[ arXiv:1209.2117 ].[175] C. de Rham,
Massive Gravity , arXiv:1401.4173 .[176] A. Einstein, Die Grundlage der allgemeinenRelativit¨atstheorie , Annalen der Physik354