Universal predictions of screened modified gravity on cluster scales
^{1} Institute of Theoretical Astrophysics, University of Oslo, Postboks 1029, 0315 Oslo, Norway
email: maxbg@astro.uio.no
^{2} Astrophysics, University of Oxford, DWB, Keble Road, Oxford, OX1 3RH, UK
Received: 27 May 2015
Accepted: 31 August 2015
Modified gravity models require a screening mechanism to be able to evade the stringent constraints from local gravity experiments and, at the same time, give rise to observable astrophysical and cosmological signatures. Such screened modified gravity models necessarily have dynamics determined by complex nonlinear equations that usually need to be solved on a modelbymodel basis to produce predictions. This makes testing them a cumbersome process. In this paper, we investigate whether there is a common signature for all the different models that is suitable to testing them on cluster scales. To do this we propose an observable related to the fifth force, which can be observationally related to the ratio of dynamicaltolensing mass of a halo, and then show that the predictions for this observable can be rescaled to a near universal form for a large class of modified gravity models. We demonstrate this using the HuSawicki f(R), the Symmetron, the nDGP, and the Dilaton models, as well as unifying parametrizations. The universal form is determined by only three quantities: a strength, a mass, and a width parameter. We also show how these parameters can be derived from a specific theory. This selfsimilarity in the predictions can hopefully be used to search for signatures of modified gravity on cluster scales in a modelindependent way.
Key words: largescale structure of Universe / dark energy / galaxies: clusters: general / galaxies: kinematics and dynamics / gravitation
© ESO, 2015
1. Introduction
A possible solution for the observed acceleration of our Universe is that general relativity (GR) is not the correct theory to describe gravity on large scales. Over the past decade, a large number of alternative gravity theories have been proposed (Clifton et al. 2012). Many of these models cannot alleviate the existing, or they introduce new problems; nevertheless, the idea that cosmic acceleration has a gravitational origin persists.
High precision tests of gravity indicate that general relativity is the correct theory for describing gravitational physics on Earth and in the solar system (Bertotti et al. 2003; Will 2006; Williams et al. 2004). This places strong constraints on any theory that seeks to modify general relativity. Thus, if gravity is modified on large scales – giving rise to cosmic acceleration – then a way of suppressing the modifications of gravity in the welltested regimes is required (Khoury 2013). A key feature of screening mechanisms currently being considered is that they are fundamentally nonlinear, thereby making the study of their cosmological effects challenging (Vainshtein 1972; Khoury & Weltman 2004a; Hinterbichler & Khoury 2010).
As previously argued in Gronke et al. (2014, 2015), deviations from GR occurring in screenedmodified gravity theories (SMGs) can generally be grouped into three categories: no deviations (fully screened regime), maximum deviations (unscreened regime), and the intermediate deviations (partially screened regime). As a natural continuation of this idea, the main questions we would like to address in this paper follow:

How do the extents and occurrences of the three regimes vary withthe screened modified gravity model and its parameters?

Is it possible to remap observables within one model and/or between different models; i.e., are SMGs selfsimilar in some way?
Previous work has been carried out in this direction. Brax et al. (2012a,b) have developed a theoretical framework in which it is possible to describe certain SMGs with two free functions. Using this framework Winther & Ferreira (2015) propose an approximate scheme to perform fast Nbody simulations of SMGs by combining linear theory with a screening factor obtained from studying spherical symmetric systems. Also, Mead et al. (2015) have recently published a technique for remapping Λ cold dark matter (CDM) power spectra to the Chameleon f(R) model with ~3% accuracy. Complementary to that, we want to base this work on an observationally focused approach. Specifically, in this paper we focus on the behavior of screened modified gravity theories inside and in the vicinity of clusters of galaxies. In general we expect the strongest signal of modified gravity to be found in cosmological voids, but clusters are easier to work with theoretically and observationally, and they remain an important observational probe of gravity on intermediate scales (Baker et al. 2015).
To answer the questions above, we have structured this paper as follows. In Sect. 2, we give a general introduction to the topography of SMGs. We also picked three models as examples, which we describe in more detail. This includes the respective solving method for halo density profiles. In Sect. 3, we present the numerical results of these three models. In addition, we study the effects of varying the model parameters on the observables and try to empirically find a rescaling that counterbalances this variation. Also in this section, we develop a general analytic rescaling method and compare it to our numerical results. In Sect. 4, we discuss our results and possible observational implications. We conclude in Sect. 5.
Throughout this work we use (Ω_{m0},H_{0}) = (0.3, 0.7), , , and denote values today with a subscript zero.
2. Methods
In this section we give a brief introduction to screened modified gravity theories. We present a few example models in more detail and then give a general framework for considering a large class of models. In the end we discuss the numerical implementation for solving for the field profiles of the models for the case of a NavarroFrenkWhite (NFW) density profile.
2.1. Modified gravity theories
When we speak about modification of gravity or modified gravity theories, we generally mean the addition of terms to the EinsteinHilbert action. Specifically, in scalartensor theories one (or multiple) additional scalar field is added with their respective Lagrangian L_{ϕ}, and the general action reads as (1)where the matter fields ψ^{(i)} couple to the metric , and R is the Ricci scalar. We limit ourselves to one scalar field and a conformal coupling to matter^{1}. Since the scalar field is coupled to matter, i.e., A = A(ϕ), an extra fifth force arises that is in the nonrelativistic limit, and per unit mass is given by (2)The predictions of general relativity have been confirmed to high accuracy in the laboratory and in the solar system, which severely constrains any new such force in our (Earth’s) local environment (Bertotti et al. 2003; Will 2006; Williams et al. 2004). If the fifth force is weak close to earth, then one would naively expect it to also be weak on cosmological scales. However, in the past decade several theories have been proposed where one is able to dynamically suppress the effect of such a fifth force in high density environments (relative to the cosmic mean) and still have interesting cosmological signatures. If this is the case, then a theory is said to possess a screening mechanism (for reviews see, e.g., Clifton et al. 2012; Khoury 2010; Joyce et al. 2015).
The action Eq. (1) allows a loose mathematical classification of possible screening mechanisms:

1.
Screening due to the scalar field value ϕ. If L_{ϕ} possesses a potential V(ϕ), it is possible to suppress the fifth force thanks to the value of the scalar field moving inside the potential. Physically, this can occur if the coupling strength depends on the local environment as is the case in the Symmetron mechanism (Hinterbichler & Khoury 2010; Hinterbichler et al. 2011). Another possibility is that the range of the fifth force depends on the environment, which is often called the Chameleon mechanism (Khoury & Weltman 2004b; Khoury & Weltman 2004a; Mota & Shaw 2007). A frequently considered SMG possessing the Chameleon screening is the Hu & Sawicki (2007)f(R) formulation.

2.
Screening due to derivatives of the scalar field ∂ϕ and/or ∂^{2}ϕ. In this case, the screening is caused by the relation between the kinetic term of L_{ϕ}, i.e., , and the rest of the scalar field Lagrangian. A popular mechanisms incorporating this idea is the Vainshtein mechanism (Vainshtein 1972) , where screening depends on the value of ∂^{2}ϕ, and the kmouflage mechanism (Babichev et al. 2009; Brax & Valageas 2014a) , where the magnitude of (∂ϕ)^{2} decides whether a region is screened or not. Particular SMGs employing the Vainshtein Mechanism are massive gravity (de Rham 2014) , Galileons (Nicolis et al. 2009) , and the DGP model (Dvali et al. 2000).
A selective overview over some models and their employed screening mechanism is given in Table 1.
Over the past decade there have been several studies of modified gravity on nonlinear scales, often using Nbody simulations. For theories that screen via the first method described above (chameleonlike screening), there have been studies of f(R) gravity (Li et al. 2012; Llinares et al. 2014; Puchwein et al. 2013; Zhao et al. 2011a; Oyaizu 2008; He et al. 2014; Lombriser et al. 2012) , other chameleon models (Brax et al. 2013) , the Symmetron (Davis et al. 2012; Llinares et al. 2014; Brax et al. 2012c) , and the environmentdependent dilation model (Brax et al. 2011, 2012c). For models that employ the second screening method (Vainshteinlike screening), there have been studies of the DGP model (Schmidt et al. 2009; Schmidt 2009; Chan & Scoccimarro 2009; Khoury & Wyman 2009; Li et al. 2013b; Falck et al. 2015, 2014) , Galileons (Barreira et al. 2013; Li et al. 2013a) , and kmoflage (Brax & Valageas 2014b). This is by no means a complete list.
In spite of the differences between the different screening mechanisms and the theories employing them, some common quantities can be defined. In particular, we focus on the enhancement of the gravitational force (3)with its theoretical maximum value γ_{max}, given a specific modified gravity model. As described in Sect. 4.2, γ can be related to the ratio of the dynamical to the lensing mass.
Overview of screened modified gravity models.
2.2. Example models
The screened modified gravity models that we considered possess three intrinsically different screening mechanisms: the Vainshtein (Vainshtein 1972) screening where the derivatives of the scalar field play a major role; the Symmetron (Hinterbichler et al. 2011a,b) screening, where the coupling strength is environment dependent and goes to zero in high density regions; and the Chameleon (Khoury & Weltman 2004a,b) screening, where the range of the fifth force approaches zero for regions with high matter density. As an example of the latter, we focus on the Hu & Sawicki (2007)f(R) model, which incorporates the Chameleon screening mechanism.
In the following, we do not present the example models in great detail. Instead, we refer the reader to our previous work (Gronke et al. 2014, 2015) , where we stated the full equations for the Symmetron and the Hu & Sawicki (2007)f(R) model or – even better – to excellent reviews, such as those by Clifton et al. (2012); Khoury (2013); Joyce et al. (2015); and, Koyama (2015).
2.2.1. Symmetron model
In the Symmetron model (Hinterbichler & Khoury 2010; Hinterbichler et al. 2011) , both the coupling function and the potential are symmetric around ϕ = 0 are given by Here, the prefactors M, μ and λ can be rewritten as which leaves the Symmetron parameters to be the scale factor at symmetry breaking a_{ssb}, the range of the fifth force in vacuum L, and the fiducial coupling β.
In a FriedmannLemaîtreRobertsonWalter (FLRW) background and in the quasistatic limit (Noller et al. 2014) , the field equation becomes (9)where . The fifth force is given by (10)In regions of high ambient matter density (), the field will reside close to φ = 0, and the fifth force will be screened. In lowdensity regions, on the other hand, the fifth force can be in full operation. In the Symmetron model, the maximum enhancement of the gravitational force, γ_{max}, is given by (11)Since β is a free parameter, the model can, in principle, give rise to an arbitrarily large (or small) deviation. How the screening mechanism works is explained in more detail in Sect. 2.3.
2.2.2. HuSawicki f(R) model
The f(R) gravity models are a group of modified gravity theories where the Ricci scalar R in the EinsteinHilbert action is replaced by a generic function R + f(R). By applying a conformal transformation, the f(R) action can be brought in the form of Eq. (1) (see, e.g., Clifton et al. 2012; Brax et al. 2008). This transformation^{2} translates f(R)gravity into a scalar tensor theory. The theory screens via the chameleon mechanism.
The fifth force is given by (12)where f_{R} ≡ df(R)/ dR is the scalar field. In a FLRW background spacetime and in the quasistatic limit (Noller et al. 2014) , the fieldequation becomes (13)where n and f_{R0} are dimensionless model parameters. The maximum enhancement of the gravitational force, i.e. the enhancement when there is no screening, is simply (14)so gravity can be enhanced by up to 33%. How the screening mechanism works is explained in more detail in Sect. 2.3.
2.2.3. DGP
Dvali et al. (2000, DGP) model is an example of a braneworld model where we are confined to live in a fourdimensional brane that itself is embedded in a fivedimensional spacetime. The DGP model can give rise to selfacceleration of the Universe without an explicit cosmological constant, but this branch has problems with instabilities. We take the normal branch (n)DGP model as our main example. To get accelerated expansion in this branch, we need to add dark energy, and for simplicity, we assume that the background expansion is the same as in ΛCDM. This choice is not expected to change any of our conclusions.
Since it is a higher dimensional theory, the DGP model only fits into the formalism of Eq. (1) as an effective theory. However, the modifications of gravity in this model can also be described as a fifth force given by (15)where the evolution of the scalar field ϕ, the cocalled branebending mode, is determined by the field equation^{3}(16)where (17)The maximum value of the gravitational force in this model is (18)For high values of r_{c}, we have β_{DGP}(a = 1) ∝ r_{c}, so The larger r_{c}, the weaker the effect of the fifth force.
As also stated in Schmidt (2010), in spherical symmetry we can integrate Eq. (16) directly to obtain (19)This translates to a solution for the fifth force profile given by (20)where (21)We have two regimes. For ϵ ≪ 1 (large r), we get and gravity is maximally modified. On the other hand, for ϵ ≫ 1 (small r), we have (22)and the fifthforce is screened.
2.3. The unified { m(a),β(a) } description of chameleonlike models
Brax et al. (2012b,a) show that any scalartensor theory that screens using the Chameleon and/or Symmetron mechanism can be described uniquely by the evolution of the mass and matter coupling along the cosmological attractor. Thus parametrizing these two functions, which have a clear intuitive meaning – the range and the strength of the fifth force, respectively – to be effectively parametrized. The mapping between the potential V(ϕ) and the coupling function A(ϕ) to { m(a),β(a) } has been derived in Brax et al. (2012b) and reads as (in parametric form) where .
One of the advantages of using this form is the direct relationship between the evolution of linear matter perturbations and the β, m functions (26)The two functions are, in the case of the Symmetron model, given by (27)i.e., they follow the same evolution. The parameters β_{0} and m_{0} are related to the model parameters described in the previous section via β_{0} = β and .
For the HuSawicki f(R) model, on the other hand, we have (28)and (29)which means that β is constant, and m(a) ∝ a^{− 3(n/ 2 + 1)} for small a. Note that the nDGP model does not map on to this parametrization except for the case of linear perturbations where it would correspond to m(a) = 0 and β(a) = 1/.
When the scalar field, in a region of density ρ is close to the minimum of the effective potential (this corresponds to the fully screened regime), we have ϕ ≈ ϕ(a_{ρ}) where . The screening condition for an object with Newtonian potential Φ_{N} in a region of density ρ can therefore be written as^{4}(30)An object is screened whenever ϵ(a_{ρ},a_{env}) ≪ 1. Here, where ρ_{env} is the density of the environment the object is located in. The min condition ensures that we get the correct value, ϵ = 1, in the nonscreened regime. The screening condition is related to γ via (31)We should note that this is a very rough analytical approximation, but it is usually good enough for orderofmagnitude estimates. To get accurate predictions, we need to solve the field equation numerically, as discussed in the next section.
Fig. 1 γ vs. R/R_{200} for NFW halos with masses ranging from 10^{9} M_{⊙}h^{1} to 10^{14.5} M_{⊙}h^{1}. The left panel shows the Symmetron solution using (a_{ssb},L,β) = (0.7, 1 Mpc h^{1}, 1), the central panel the Hu & Sawickif(R) solution with ( f_{R0}  ,n) = (10^{6}, 1), and the right panel the nDGP solution with r_{c} = 1. 

Open with DEXTER 
2.4. Numerical solving methods
To obtain accurate predictions for the field and force profiles in the models we study we need to numerically solve the field equations Eqs. (9) and (13). We solve these equations by discretizing them on a grid, using the NFW density profile for ρ_{m} (see Appendix A for the relevant NFWequations), and then use NewtonGauss Seidel relaxation with multigrid acceleration to obtain the solution.
In spherical symmetry the field equation becomes onedimensional, φ = φ(r) and . For the Symmetron model, we implemented a simple fixedgrid multigrid solver for this onedimensional problem. To do this we write Eq. (9) as ℒ_{i} = 0 where (32)where subscript i denotes the value at the ith gridpoint and (33)is a secondorder discretization of the Laplacian. We start by making a guess for the solution, and then we loop through the grid updating the solution using (34)The multigrid technique is used to speed up convergence. For more details about the methods we have used see Llinares et al. (2014).
Fig. 2 Semianalytical solutions for γ for the Symmetron (left) and f(R) (right) using the same parameter values (and color coding) as in Fig. 1. 

Open with DEXTER 
The f(R) equation Eq. (13) is stiffer than the Symmetron equation for low values of  f_{R0} , and a simple fixedgrid onedimensional solver does not converge for the whole range of considered parameters. The main reason for this is that  f_{R}  is required to be strictly larger than zero (but can get very close). This means that one has to ensure that  f_{R}  does not become negative in the solving process. To achieve this, we used the already welltested f(R)module in the ISIS code (Llinares et al. 2014). The only modification to the module is the change to a purely geometrical refinement criterion based on the value of ρ_{m}. In total, we used 10−15 levels of refinements starting from a base grid of N = 64^{3} grid nodes in order to get accurate field profiles down to very small radii. The Symmetron is also implemented in the ISIS code, and we tested that our Symmetron field profiles agree very well with those obtained with our code.
After having calculated the field profiles, the fifth force can be calculated by using Eqs. (10) and (12) for the Symmetron and f(R) model, respectively. For nDGP the analytical solutions to the field equation and force profile in spherical symmetry is given by Eqs. (19) and (20) , respectively. This makes it needless to solve any differential equations numerically^{5}. The derivation of the analytical solution is shown in Appendix B.
3. Results
In this section we present our results from numerical solutions of the field equations (Sects. 3.1, 3.2) and the obtained rescaling method for the three example models. The rescaling is first found empirically in Sect. 3.3 and then semianalytically in the { m(a),β(a) } formulations in Sect. 3.4.
Fig. 3 vs halo mass M. The left panel shows the Symmetron curves with (a_{ssb},L,β) being (from solid to more interceptions) 0.7, 1 Mpc h^{1}, 1 (red), 0.7, 2 Mpc h^{1}, 1 (blue), 0.3, 1 Mpc h^{1}, 1 (green), 0.7, 1 Mpc h^{1}, 2 (purple). The central panel shows the Hu & Sawickif(R) solution with log _{10}  f_{R0}  = (−4, −5, −6, −7) (in red, blue, green, and purple, respectively) and the right panel the nDGP solution with r_{c} = (0.5, 1 ,3) (from top to bottom). 

Open with DEXTER 
3.1. Force profiles
Figure 1 shows the radial profiles of γ obtained solving Eqs. ((9) , (13) , (20) ) for the Symmetron, f(R), and, nDGP model, respectively (from left to right panels). We display a set of halo masses for each model using fixed model parameters.
The differences in shape between the models are striking. The Symmetron model shows a whole range of solutions: from fully screened in the center for the heaviest clusters to not screened at all for the lightest ones. The maximum of the fifth force seems to be located ~1−10R_{200} and drops again for larger radii. This drop happens due to the fifth force having a finite range, .
The f(R) model presented in the central panel of Fig. 1 inherits a much more abrupt change from fully screened to not screened in the central halo region. With this particular choice of parameters ( f_{R0}  = 10^{6},n = 1), halos with a mass of ≳10^{12} M_{⊙} seem to be fully screened in the central (R ≲ 0.1R_{200}) region, whereas in lighter halos, the fifth force stays at its theoretical maximum 1/3F_{N}. For comparison, in Fig. 2 we show the force profiles derived using the semianalytical method described in Sect. 3.4 for the f(R) (left) and Symmetron (right panel) models. The agreement is fairly good around the most significant region R ~ 1−5R_{200} where γ peaks. The semianalytic results are presented in more detail in Sect. 3.4.
In the righthand panel of Fig. 1, we show the radial scalar field profile ϕ(r) for a NFW halo in the nDGP model with r_{c} = 1. As is also possible to see from Eq. (B.2) , the solution includes a curious feature: owing to the concentrationmass relation used (Eq. (A.3) ), heavier halos are less screened in the central region of the halo. Generally, the fifth force is less strong in the nDGP model and follows a similar functional form that is almost independent of mass^{6}.
3.2. Profile variation with halo mass
To capture the change of the fifth force profile with halo mass, we chose to show the massweighted mean of the fifth force (35)where w(x) is the normalized mass fraction w(x) dx ≡ dM/M( <x_{cut}) as defined in Eq. (A.4). This quantity was introduced by Schmidt (2010) for the f(R) and DGP model and also considered in Clampitt et al. (2012) for the symmetron model. The outer cutoff was fixed to x_{cut} = 10. This value was chosen so that w(x_{cut}) is roughly onetenth of the maximum weight, meaning that not too much information is lost. As discussed in Sect. 4.2 (see Eq. (55) ), can be related to the ratio of the dynamical mass to the lensing mass.
Figure 3 shows versus the halo mass. However, this measure is not unique. Instead it is also possible, for example, to use a difference between the force in the center and in the outskirts of the halo. We discuss the advantages and caveats of using as an “observable” in Sect. 4.2.
The lefthand panel of Fig. 3 shows three Symmetron solutions with different model parameters (a_{ssb},L,β). The slope, as well as the offset, of the three curves is seen to be fairly different. Also notable is the immense change in with changing halo mass. The f(R) solutions (shown in the central panel of Fig. 3) possess an upper limit at the already mentioned value of onethird. For lighter halos falls off rapidly. In the rightermost panel of Fig. 3, we show the variation in with mass for the nDGP model for three different values of r_{c}. Overall, the impact of the fifth force is much weaker than in the other models. As discussed in Sect. 3.1, the variation in the fifth force with halo mass is only due to the massconcentration relation of the NFW halo. This feature manifests itself in Fig. 3 in the slight increase in with increasing halo mass.
3.3. Empirical rescaling
To find a rescaling that maps the various characteristics presented in the last section onto one single curve, we need to rescale the x and the y axes of Fig. 3. The latter is naturally normalized to the theoretical possible maximum of the gravitational enhancement γ_{max}. The rescaling of the halo mass is less trivial. Therefore, we simply define a characteristic mass μ_{200} in the following way: (36)where is still the massweighted average of the enhancement of the fifth force as defined in Eq. (35) In general μ_{200} is a function of the various model parameters. Halos with mass much higher than μ_{200} are in the fully screened regime, whereas halos with mass much lower than μ_{200} are in the unscreened regime.
In this section we find the functional form of μ_{200} for each model merely empirically from our numerical solutions. However, in the next section, we derive semianalytic predictions for a large class of scalar tensor theories.
This means that, by looking at Fig. 3, it is possible to read off several values of μ_{200}. For example, for the f(R) model with (n,  f_{R0}  ) = (1, 10^{6}) (in green), γ( ~10^{13} M_{⊙} h^{1}) = 1/6 = 1/2γ_{max}, i.e., μ_{200}, is in that particular case ~10^{13} M_{⊙} h^{1}.
For the Symmetron model, we found μ_{200} to be well fit by (37)and γ_{max} is given by Eq. (11). Figure 4 shows the resulting rescaled curves, i.e., versus M/μ_{200}. In addition to the four parameter sets already presented in the lefthand panel of Fig. 3 (keeping the color coding), we plot an additional 125 curves with randomly chosen parameters. Specifically, we drew L ~ [0.2, 2] Mpc h^{1}, a_{ssb} ~ [0.2, 0.8], and β ~ [0.2, 2.4] (all uniformly).
It is notable that in spite of this wide parameter range, the resulting curves in Fig. 4 are very similar. All follow the characteristic shape – from γ ~ γ_{max} at M ≪ μ_{200} through γ(M ≈ μ_{200}) = γ_{max}/ 2 to γ ~ 0 for M ≫ μ_{200} – with very little scatter. We note, however, that for low masses, does not fully approach γ_{max} but falls ~10% short. This is due to the drop in γ for larger radii as is visible in the lefthand panel of Fig. 1.
Fig. 4 versus M/μ_{200} for the Symmetron model. The solid colored lines indicate the same model parameters as in Fig. 3. The gray lines are an additional 124 randomly chosen Symmetron parameters. μ_{200} for the Symmetron is given by Eq. (37). The black dashed lines indicate the point defining μ_{200} in Eq. (36) , and the gray arrows point to the screened and unscreened regimes. 

Open with DEXTER 
Fig. 5 versus M/μ_{200} for the HuSawicki f(R) model. The solid colored lines indicate the same model parameters as in Fig. 3. The gray lines are an additional 5 values of f_{R0}. μ_{200} for this model is given by Eq. (38). 

Open with DEXTER 
In Fig. 5 we present the results of the same rescaling technique applied to the f(R) results. In this case the rescaling mass is found to be (38)and γ_{max} = 1/3. The curves in Fig. 5 correspond to the numerical solutions for each combination of n = (1, 2) and  f_{R0}  = (10^{8}, 3 × 10^{8}, 10^{7}, 3 × 10^{7},...,10^{4}), with the ones already presented in Fig. 3 colorcoded accordingly. Also here, the scatter around the mean form is quite small. The second striking feature the similarity in the overall shape of the curves in Figs. 4 and 5 – even though the solution stems from two independent modified gravity theories. The main differences between the two models are the slightly different transitions widths and, the fact that for M ≪ μ_{200}, γ ≈ γ_{max} in the f(R) case, whereas the Symmetron curves are limited to a lower value.
Rescaling the mass relation for the nDGP model (presented in the right panel of Fig. 3) is not as straightforward. Since the curves are nearly constant (and do not possess a characteristic drop), no empirical formulation for μ_{200} is apparent. However, since calculating in the nDGP model is not computationally expensive, we can numerically solve Eq. (36) for several values of r_{c}. This let us fit the empirically found relationship (39)which yields (A,a,b) = (11.14, 6.72, −0.57). The rescaled solutions obtained with r_{c} = (0.5, 1,..., 20) showed so little variation that we chose not to explicitly show them. Instead, we display in Fig. 6 the average curves (solid lines), as well as the maximal deviation found for each (dashed lines). Clearly, the nDGP curves show very little variation.
Fig. 6 versus M/μ_{200} for the three example models. The solid [dashed] lines show the mean (minimum and maximum) values for the 128 (Symmetron), 16 (f(R)), and 20 (nDGP) numerical curves obtained. The solid gray lines show the fits of Eq. (41) described in Sect. 3.3. 

Open with DEXTER 
Figure 6 once more shows the difference between the Symmetron and f(R) model with their characteristic “inverseSshape” on the one side, and, the nDGP – with a loglinear relation – on the other. A simple fitting function that is able to match the characteristic “inverseSshape” given by (40)However, it is also clear that the models of the former group do not show exactly identical behavior. Instead, the transition region for the Symmetron model is a bit wider, and, as mentioned before, it does not converge to γ_{max} for lower masses. This is why we propose the more general form (41)with two free parameters , α. Here, the prefactor can be interpreted as the part of the theoretical enhancement of gravity that is effectively active (and, thus, can only be <1), and α controls the width of the partially screened region with α = 0 being “infinitely wide” and α → ∞ the transition turning into a step function. Also, if α< 0, the “inverseSshape” flips vertically and turns into an “Sshape”. Quantitatively, we can relate α to the half width of the partially screened region through (42)Here, M_{200} ∈ [μ_{200}/W,Wμ_{200}] is the interval where .
Figure 6 also shows the results of fitting Eq. (41) to all of our numerical data (with enforcing ). It is noteworthy that this simple, twoparameter fit describes the numerical results of the three studied models very well, even for the nDGP model. As expected, the α parameter reflects the fundamental difference between the nDGP and the two other models, as it is close to zero (−0.03) for the nDGP and more similar for the other two models. Specifically, we obtained α = 0.46 and 0.78 for the Symmetron and f(R) model, respectively.
3.4. Semianalytical derivation of μ_{200} in the { m(a),β(a) } formulation
In the following, we discuss the origin of the empirical rescalings found. Using the general framework of the { β(a),m(a) } formulation described in Sect. 2.3, we can derive an analytic approximation for μ_{200}. It turns out that we can do this quite generally, so our results below can be applied to any scalartensor theory with a Lagrangian and a conformal coupling to matter. We only focus on a very rough derivation, which can be seen as the foundation of a more precise treatment in the future.
Before we do this, we should check to what extent our semianalytic approximation, which formally holds only for a spherical tophat overdensity, is valid for NFW halos. In Fig. 2 we show the semianalytical force profiles that can be compared to Fig. 1, which shows the true numerical results. The qualitative agreement is fairly good. We are able to match the shape around the peak, which is the significant part of the profile that will give rise to most of the signal in , and the amplitude is off by no more than a factor of a few, which is good enough for our purposes.
As a first – rather crude – approximation, we assume that a halo is fully screened within a certain radius x_{screen}. In the unscreened regime x_{screen} → 0, but generally 0 ≲ x_{screen} ≲ 10. Thus, we approximate γ(x) ≃ γ_{max}θ(x_{screen}−x) where θ is the Heaviside function. This leads to a solution for the massweighted average of γ, which reads as (43)Furthermore, using the definition of μ_{200}, namely , we require (44)which requires for typical halo masses.
The dependence of x_{screen}, hence of μ_{200}, on the various SMG model parameters can be calculated using the screening condition Eq. (30). In this case, the condition reads as With and M_{200} = 4π/ 3Δ_{vir}Ω_{m}ρ_{c0}r^{3} → r ∝ M^{1/3}, we can easily calculate how μ_{200} scales with the model parameters in any given model.
For the Symmetron we find the scaling M( <. As a result, μ_{200} scales as (47)where the proportionality factor was calculated numerically. This scaling is exactly the same as we found empirically in the previous section, but the proportionality factor is off by a factor of ~3−10.
Fig. 7 Semianalytical predictions for versus M/μ_{200} for the HuSawicki f(R) model (red), the Symmetron (blue) and the Dilaton model (purple). The solid black curve shows the fitting formula Eq. (40). The shaded regions show the results of the full numerical analysis. 

Open with DEXTER 
For the HuSawicki f(R) model, we find the scaling μ_{200} ∝  f_{R0}  ^{3/2}C^{3/2(n + 1)} M_{⊙}/h, where C is a constant of order one, dependent only on x_{screen} and Ω_{m}. Numerically we find that the results are only weakly dependent on n: changing n from 1 to 5 only gives us a factor 1.5 offset in the value we find for μ_{200}. Again the proportionality factor is found numerically and gives us (48)which is almost spot on the empirical relation found in the previous section.
For a general { m(a),β(a) } model, we see from Eq. (23) that we would (again very roughly) expect to have . This in the condition Eq. (45) means that we generally expect where λ_{φ}(a) = m^{1}(a) is the range of the fifth force.
In Fig. 7 we show the semianalytical profiles for in terms of M_{200}/μ_{200} for f(R) and the symmetron, together with the full numerical results. To demonstrate that this rescaling works in general, and not just for the two models presented here, we also calculated the predictions for the Dilaton model presented in Brax et al. (2012c, 2010). This model is characterized by where r,s,m_{0}, and β_{0} are free model parameters and where the evolution of β are very different from our two example models. Figure 7 shows that the rescaling also works perfectly for the Dilaton model.
We have not been able to get a prediction for the width of the partially screened region semianalytically, so we leave this for future work.
Fig. 8 γ(M) for NFW halos with masses ranging from 10^{9} M_{⊙}h^{1} to 10^{15} M_{⊙}h^{1} using Eq. (52) with . The four bands of lines show the results (from top to bottom) of ζ = 0.1,1.0,10,100, and within each band, we have (from top to bottom) n = 0.5,0.4,0.3, and 0.2. 

Open with DEXTER 
3.5. Predictions for other screening mechanisms
For other screening mechanisms that are not encapsulated by the { m(a),β(a) } formulation, it is harder to make general predictions. For any theory with a derivative shift symmetry (the Galileon symmetry), secondorder equations of motion the field equation in spherical symmetry are integrable, leading to an algebraic equation on the form (51)where the last fraction can be identified as g_{NFW}(cx) /g_{NFW}(x) (see Eq. (A.2) ). For a general model this leads to (52)for some function G that satisfies G(0) = 1 and G(∞) = 0. For DGP we have with .
For models that only have a shift symmetry, like the kmouflage screening mechanism, the field equation, (53)can also be integrated up, and from this it follows that the fifthforce F_{φ} ∝ ∇φ is a function of for spherical symmetry. This leads to a slightly different form (54)for some function F satisfying F(∞) = 0.
The exact form for the functions F and G are model dependent, making it hard to make definite statements here. However, one general feature we are able to deduce is that, just as we saw for DGP, γ is only sensitive to the concentration of a halo. We therefore expect a fairly weak mass dependence of for models in this class.
As an example, we can try to parametrize the functional form for F and G and calculate . One functional form that satisfies the requirements is with n ∈ (0,1). The results we get from this of course depend on the values we take for the parameters (like ζ); however, taking realistic values for the parameters, which is motivated by the requirement of satisfying local gravity constraints, the typical result we find is as expected (see Fig. 8) : a very weak M dependence of .
Fig. 9 Rescaling parameter μ_{200} versus the maximum gravitational enhancement γ_{max} for several discussed models (solid colored lines) and the respective partially screened regions (colored semitransparent areas). The red line shows the space spanned by the HuSawicki f(R) model when varying  f_{R0}  from 10^{8} (left end) to 10^{4} (right end), while leaving n = 1 fixed. For the Symmetron, we also varied one parameter per curve (in blue, green, and violet), while leaving the others fixed at β = 1, L = 1 Mpc h^{1}, a_{ssb} = 1/2. The orange line indicates the region for nDGP models with r_{c} ∈ [0.5, 10] (left to right) with the associated semitransparent area spanning wider than the plot range. See Sect. 4.1 for details. 

Open with DEXTER 
4. Discussion
Our numerical force profiles presented in Sects. 3.1 and 3.2 are in good agreement to what has been previously found in Nbody simulations (e.g., Zhao et al. 2011a; Falck et al. 2015). The simple rescaling procedure based on the maximum enhancement of the gravitational force γ_{max} and the intermediately screened mass μ_{200} (Sect. 3.3) demonstrates that in spite of their foundational differences, the three screenedmodified gravity models investigated show selfsimilar behavior in cluster environments.
4.1. Observational implications
Our rescaling methods allow an abstract characterization of screened modified gravity theories, which themselves might possess quite different screening mechanisms. This potentially has interesting implications for constraining SMGs.
The colored lines in Fig. 9 show γ_{max} versus μ_{200} for some selected models. The lines were drawn by altering one model parameter while leaving the others constant (see figure caption for details). Each set of SMG model parameters occupies a specific point on this graph; for example, the combination (a_{ssb},β,L) = (1/2, 1, 1 Mpc h^{1}) maps to (γ_{max},μ_{200}) ≈ (1.75, 4.5 × 10^{11} M_{⊙}h^{1}). The semitransparent areas in Fig. 9 denote the extent of the partially screened region using Eq. (42) with d = 0.2. This means that for every selected value of μ_{200}, the semitransparent area indicates the range where the fifth force is at least 20% but a maximum of 80% active.
Thus, Fig. 9 indicates (i) what halos are partially screened and; (ii) what the maximum enhancement of gravity is for this model. Physically, μ_{200}/W can be interpreted as the minimum mass where one expects a variation in halo properties; that is to say:

For M_{200} ≲ μ_{200}/W halos in this model are unscreened, thus allowing for comparison with the ΛCDM predictions. In particular, one expects that the majority of halo properties are rescaled by ~γ_{max}. This can make the comparison with ΛCDM difficult when using, for example, the dynamical mass as mass estimate.

Halos with M_{200} ∈ [μ_{200}/W,Wμ_{200}] are partially screened. This means that in this mass range, halo properties might differ from the expected. Since halos in this mass range are between screened and unscreened, the environment of these halos has a major effect (see Zhao et al. 2011b; Winther et al. 2012, for the study of this effect using Nbody simulations). Also, comparing the inner and outer regions of individual objects can lead to insight (for instance, the ratio M_{500}/M_{200}). In Fig. 9 the partially screened region is shown as a shaded area in the matching color. As demonstrated for the nDGP case, it is possible that – for some SMGs – all halos are always partially screened.

For masses ≳Wμ_{200}, the halo is fully screened, so no deviation from the ΛCDM predictions is expected.
Conclusively, in any observablemass scaling relation, we predict a peculiar feature at μ_{200}, where the observables match the ΛCDM prediction for higher masses and where they are different for lower masses. However, the latter statement only holds if the true mass of the halo is used. This is because mass tracers like the dynamical mass can also be affected by the fifth force (Schmidt 2010).
The extent of the partially screened region, which we quantified with α (or, W, see Eq. (42) , can potentially be used to distinguish between various screening mechanisms. We found that this extent is constant for each screenedmodified gravity considered, but differs significantly between them. In particular, the value of α found in the numerical fits of Sect. 3.3 translates to log _{10}W ~ 16, ~1, and ~0.6 for the nDGP, Symmetron, and the f(R) models, respectively (using d = 0.2, as above). In a γ_{max} α  (or equivalently, a γ_{max}W) diagram, we expect eventual constraints to rule out the area γ_{max}>γ_{constraint},  α  ≲ α_{constraint} (W ≳ W(α_{constraint})).
4.2. Caveats
Several assumptions went into deriving the presented formalism. Here, we want to mention and investigate them.

Halo profiles and shape. Throughout this work, we assumed spherical symmetry and perfect NFW halos.However, in reality clusters of galaxies do not exactlyfollow an NFW mass distribution, and halos also possess a nonzero ellipticity (e.g., Oguriet al. 2010). Another phenomenonfalling in the same category is our assignment of a unique massconcentration mapping, whereas the spread in this relation isfairly wide (Neto et al. 2007). Bothare consequences of the fact that clusters are individual objectswith their own histories and environments. As a result, individual clusters (as any astrophysical object) have mostly a limitedexplanatory power and must be stacked for analysis.

Model assumptions. To rescale the curves accordingly, an intermediatemass scale has to exist where . Since this is the scale between the screened and unscreened regimes and, our work addresses “screened modified gravity theories”, we think that this can be safely assumed. The universal rescaling works perfectly for the models we have tested, and we conjecture that it will hold for all { m(a),β(a) } models. However, this might not be true for cases where m(a) and/or β(a) are very steep functions of a, leading to a very rapid evolution of the screening condition with density.

Choice of parameters. Throughout this work, we fixed a number of parameters. Probably, the most influential ones are the halo cutoff radius x_{cut} = 10 and defined μ_{200} as with k = 1/2. The latter choice is arguably quite natural. In spite of that, we do not expect our results to change dramatically if a (slightly) higher or lower value of k is taken since the model curves are expected to shift equally. The choice of x_{cut} affects our results more drastically, as can be inferred from Fig. 1. A lower value of x_{cut} leads, for instance, to a decrease in the extent of the partially screened region (i.e., an increase in α). Nevertheless, since all models are affected by this, the precise choice of x_{cut} (within reason) does not alter our conclusions.

Practicality. We focused merely on the force profiles directly, leaving the question how practical the study is since γ(r) and are not directly observable. Although this is not wrong, we caution that the fifth force triggers many observable phenomena, as has been shown in various Nbody simulations. In particular, γ(r) can be related to the ratio of the dynamical mass to the lensing mass (Schmidt 2010; Zhao et al. 2011b; Winther et al. 2012) (55)where Here, (Φ + Ψ)/2 is the dynamical potential, (Φ−Ψ)/2 is the lensing potential, and Φ and Ψ are the two metric potentials. For the models considered in this paper, the lensing potential is the same as in ΛCDM, but there are several models where the lensing potential is modified as in the Galileon, for example^{7}. In practice, however, it is not possible to measure the potentials directly. Instead, the dynamical mass is usually obtained via the cluster velocity dispersion through a powerlaw relation obtained from (ΛCDM) Nbody simulations (e.g., Becker et al. 2007). The lensing mass can be read from the stacked density profiles assuming a constant massrichness relation (Johnston et al. 2007; Hoekstra et al. 2013). This means that our massweighted average of the gravitational enhancement is directly proportional to the ratio of the dynamical to the lensing mass as has been shown by Schmidt (2010)^{8}.
In conclusion, our work can be seen as a first step toward a truly unified description of SMGs on cluster scales. Naturally, this first step can be tuned further and extended to improve its precision and include more models.
5. Conclusions
Starting from the simple idea that screened modified gravity theories in general possess three regimes within nonlinear astrophysical structures, a fully screened, an unscreened, and a partially screened regime, we investigated whether one can specify, independently of the SMG model, at which halo mass range these regimes are located and how extended they are.
These two questions can be answered as follows:

It is possible to formulate the mass of the partially screened regimeas a function of various model parameters. We demonstrated thisempirically for the Symmetron, the HuSawicki f(R), and the nDGP models (Sect. 3.3). In addition, we found a semianalytic derivation of this functional form in the { m(a),β(a) }parametrization of Brax et al. (2012b). This allowed us to expand our method easily to the Dilaton model.

The extent of the partially screened regime is an intrinsic property of a particular screening mechanism. We found that the mass range of the partially screened region differs significantly between the three example models but are (approximately) constant within a particular model.
In spite of the differences in the nature of their screening mechanisms (see Sect. 2.1) , we characterize the Symmetron, the HuSawicki f(R), and the nDGP models with three parameters: (i) the center of the partially screened regime μ_{200}; (ii) the maximal enhancement of the fifth force γ_{max}; and the (iii) effectiveness of the screening α (which controls the extent of the partially screened regime). We derived the first two parameters analytically in the { m(a),β(a) }framework, which includes the Symmetron, the HuSawicki f(R), and the Dilaton models.
Quantifying observational constraints in the proposed (μ_{200},γ_{max},α)parametrization thus allows several screened modified gravity models to be constrained in bulk instead of having to consider them – and their respective parameter space – individually.
Disformal coupling allows for additional screening mechanisms (Koivisto et al. 2012).
While this paper was in the final stages of completion, Barreira et al. (2015) showed that for the Cubic Galileon model, the lensing potential within clusters is practically unmodified.
In fact, Schmidt (2010) derive the relation (see their Eq. (64).).
Acknowledgments
We thank the anonymous reviewer for useful comments. M.G. thanks the Astrophysics department at the University of Oxford for their hospitality. D.F.M. acknowledges support from the Research Council of Norway through grant 216756. H.A.W. is supported by the BIPAC and the Oxford Martin School.
References
 Babichev, E., Deffayet, C., & Ziour, R. 2009, Int. J. Mod. Phys. D, 18, 2147 [NASA ADS] [CrossRef] [Google Scholar]
 Baker, T., Psaltis, D., & Skordis, C. 2015, ApJ, 802, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Barreira, A., Li, B., Hellwing, W. A., Baugh, C. M., & Pascoli, S. 2013, J. Cosmol. Astropart. Phys., 10, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Barreira, A., Li, B., Jennings, E., et al. 2015, arXiv eprints [arXiv:1505.03468] [Google Scholar]
 Becker, M. R., McKay, T. A., Koester, B., et al. 2007, ApJ, 669, 905 [NASA ADS] [CrossRef] [Google Scholar]
 Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Brax, P., & Valageas, P. 2014a, Phys. Rev. D, 90, 023507 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., & Valageas, P. 2014b, Phys. Rev. D, 90, 123521 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., van de Bruck, C., Davis, A.C., & Shaw, D. J. 2008, Phys. Rev. D, 78, 104021 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Brax, P., van de Bruck, C., Davis, A.C., & Shaw, D. 2010, Phys. Rev. D, 82, 063519 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., van de Bruck, C., Davis, A.C., Li, B., & Shaw, D. J. 2011, Phys. Rev. D, 83, 104026 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Davis, A.C., & Li, B. 2012a, Phys. Lett. B, 715, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Davis, A.C., Li, B., & Winther, H. A. 2012b, Phys. Rev. D, 86, 044015 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Davis, A.C., Li, B., Winther, H. A., & Zhao, G.B. 2012c, J. Cosmol. Astropart. Phys., 10, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Davis, A.C., Li, B., Winther, H. A., & Zhao, G.B. 2013, J. Cosmol. Astropart. Phys., 4, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. C., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 104005 [NASA ADS] [CrossRef] [Google Scholar]
 Clampitt, J., Jain, B., & Khoury, J. 2012, J. Cosmol. Astropart. Phys., 1, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Davis, A.C., Li, B., Mota, D. F., & Winther, H. A. 2012, ApJ, 748, 61 [NASA ADS] [CrossRef] [Google Scholar]
 de Rham, C. 2014, Liv. Rev. Relat., 17, 7 [Google Scholar]
 Dvali, G., Gabadadze, G., & Porrati, M. 2000, Phys. Lett. B, 485, 208 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Falck, B., Koyama, K., Zhao, G.b., & Li, B. 2014, J. Cosmol. Astropart. Phys., 7, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Falck, B., Koyama, K., & Zhao, G.B. 2015, J. Cosmol. Astropart. Phys., 7, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Gronke, M. B., Llinares, C., & Mota, D. F. 2014, A&A, 562, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gronke, M., Llinares, C., Mota, D. F., & Winther, H. A. 2015, MNRAS, 449, 2837 [NASA ADS] [CrossRef] [Google Scholar]
 He, J.h., Li, B., Hawken, A. J., & Granett, B. R. 2014, Phys. Rev. D, 90, 103505 [NASA ADS] [CrossRef] [Google Scholar]
 Hinterbichler, K., & Khoury, J. 2010, Phys. Rev. Lett., 104, 231301 [NASA ADS] [CrossRef] [Google Scholar]
 Hinterbichler, K., Khoury, J., Levy, A., & Matas, A. 2011, Phys. Rev. D, 84, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Hoekstra, H., Bartelmann, M., Dahle, H., et al. 2013, Space Sci. Rev., 177, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 1 [Google Scholar]
 Johnston, D. E., Sheldon, E. S., Wechsler, R. H., et al. 2007, arXiv eprints [arXiv:0709.1159] [Google Scholar]
 Joyce, A., Jain, B., Khoury, J., & Trodden, M. 2015, Phys. Rep. 568, 1 [Google Scholar]
 Khoury, J. 2010, arXiv eprints [arXiv:arXiv:1011.5909] [Google Scholar]
 Khoury, J. 2013, arXiv eprints [arXiv:arXiv:1312.2006] [Google Scholar]
 Khoury, J., & Weltman, A. 2004a, Phys. Rev. D, 69, 044026 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Khoury, J., & Weltman, A. 2004b, Phys. Rev. Lett., 93, 171104 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Khoury, J., & Wyman, M. 2009, Phys. Rev. D, 80, 064023 [NASA ADS] [CrossRef] [Google Scholar]
 Koivisto, T. S., Mota, D. F., & Zumalacárregui, M. 2012, Phys. Rev. Lett., 109, 241102 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K. 2015, arXiv eprints [arXiv:1504.04623] [Google Scholar]
 Li, B., Zhao, G.B., Teyssier, R., & Koyama, K. 2012, J. Cosmol. Astropart. Phys., 1, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Barreira, A., Baugh, C. M., et al. 2013a, J. Cosmol. Astropart. Phys., 11, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Zhao, G.B., & Koyama, K. 2013b, J. Cosmol. Astropart. Phys., 5, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., Mota, D. F., & Winther, H. A. 2014, A&A, 562, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lokas, E. L., & Mamon, G. A. 2001, MNRAS, 321, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Lombriser, L., Koyama, K., Zhao, G.B., & Li, B. 2012, Phys. Rev. D, 85, 124054 [NASA ADS] [CrossRef] [Google Scholar]
 Mead, A. J., Peacock, J. A., Lombriser, L., & Li, B. 2015, MNRAS, 452, 4203 [NASA ADS] [CrossRef] [Google Scholar]
 Mota, D. F., & Shaw, D. J. 2007, Phys. Rev. D, 75, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1995, MNRAS, 275, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 1, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450 [NASA ADS] [CrossRef] [Google Scholar]
 Nicolis, A., Rattazzi, R., & Trincherini, E. 2009, Phys. Rev. D, 79, 064036 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Noller, J., von BraunBates, F., & Ferreira, P. G. 2014, Phys. Rev. D, 89, 023521 [NASA ADS] [CrossRef] [Google Scholar]
 Oguri, M., Takada, M., Okabe, N., & Smith, G. P. 2010, MNRAS, 405, 2215 [NASA ADS] [Google Scholar]
 Oyaizu, H. 2008, Phys. Rev. D, 78, 123523 [NASA ADS] [CrossRef] [Google Scholar]
 Puchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F. 2009, Phys. Rev. D, 80, 123003 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F. 2010, Phys. Rev. D, 81, 103002 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F., Vikhlinin, A., & Hu, W. 2009, Phys. Rev. D, 80, 083505 [NASA ADS] [CrossRef] [Google Scholar]
 Vainshtein, A. I. 1972, Phys. Lett. B, 39, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Will, C. M. 2006, Liv. Rev. Relat., 9, 3 [Google Scholar]
 Williams, J. G., Turyshev, S. G., & Boggs, D. H. 2004, Phys. Rev. Lett., 93, 261101 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., Mota, D. F., & Li, B. 2012, ApJ, 756, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., & Ferreira, P. G. 2015, Phys. Rev. D, 91, 123507 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, G.B., Li, B., & Koyama, K. 2011a, Phys. Rev. D, 83, 044007 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, G.B., Li, B., & Koyama, K. 2011b, Phys. Rev. Lett., 107, 071303 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: NFW equations
The NavarroFrenkWhite (NFW) profile is an empirically derived halo profile (Navarro et al. 1995, 1996; 1997) from Nbody simulations. The NFW density profile reads as a function of normalized distance from the center x ≡ r/R_{vir}: (A.1)where (A.2)Here, c is the concentration parameter, which we link to M_{200} through the observed relationship found in Bullock et al. (2001): (A.3)Throughout this work we define Δ_{vir} ≡ 200, which means R_{vir} = R_{200} and .
Although the density of the NFW profile is diverging for x → 0, the enclosed mass within a certain radius does not (Lokas & Mamon 2001). This leaves the mass fraction within [x ± dx/ 2] as (A.4)where x_{cut} is defined below Eq. (35). This will be useful when considering massweighted quantities.
To calculate the screening condition we also need the expression for Φ_{N}. For a NFW halo, we find (A.5)
Appendix B: DGP solutions for NFW halos
In this appendix we present the equations for γ for the DGP model.
For a NFW density profile (Eq. (A.1) ) halo, we have (B.1)where . This leads to the fifth force profile being given by (B.2)where and (B.3)
All Tables
All Figures
Fig. 1 γ vs. R/R_{200} for NFW halos with masses ranging from 10^{9} M_{⊙}h^{1} to 10^{14.5} M_{⊙}h^{1}. The left panel shows the Symmetron solution using (a_{ssb},L,β) = (0.7, 1 Mpc h^{1}, 1), the central panel the Hu & Sawickif(R) solution with ( f_{R0}  ,n) = (10^{6}, 1), and the right panel the nDGP solution with r_{c} = 1. 

Open with DEXTER  
In the text 
Fig. 2 Semianalytical solutions for γ for the Symmetron (left) and f(R) (right) using the same parameter values (and color coding) as in Fig. 1. 

Open with DEXTER  
In the text 
Fig. 3 vs halo mass M. The left panel shows the Symmetron curves with (a_{ssb},L,β) being (from solid to more interceptions) 0.7, 1 Mpc h^{1}, 1 (red), 0.7, 2 Mpc h^{1}, 1 (blue), 0.3, 1 Mpc h^{1}, 1 (green), 0.7, 1 Mpc h^{1}, 2 (purple). The central panel shows the Hu & Sawickif(R) solution with log _{10}  f_{R0}  = (−4, −5, −6, −7) (in red, blue, green, and purple, respectively) and the right panel the nDGP solution with r_{c} = (0.5, 1 ,3) (from top to bottom). 

Open with DEXTER  
In the text 
Fig. 4 versus M/μ_{200} for the Symmetron model. The solid colored lines indicate the same model parameters as in Fig. 3. The gray lines are an additional 124 randomly chosen Symmetron parameters. μ_{200} for the Symmetron is given by Eq. (37). The black dashed lines indicate the point defining μ_{200} in Eq. (36) , and the gray arrows point to the screened and unscreened regimes. 

Open with DEXTER  
In the text 
Fig. 5 versus M/μ_{200} for the HuSawicki f(R) model. The solid colored lines indicate the same model parameters as in Fig. 3. The gray lines are an additional 5 values of f_{R0}. μ_{200} for this model is given by Eq. (38). 

Open with DEXTER  
In the text 
Fig. 6 versus M/μ_{200} for the three example models. The solid [dashed] lines show the mean (minimum and maximum) values for the 128 (Symmetron), 16 (f(R)), and 20 (nDGP) numerical curves obtained. The solid gray lines show the fits of Eq. (41) described in Sect. 3.3. 

Open with DEXTER  
In the text 
Fig. 7 Semianalytical predictions for versus M/μ_{200} for the HuSawicki f(R) model (red), the Symmetron (blue) and the Dilaton model (purple). The solid black curve shows the fitting formula Eq. (40). The shaded regions show the results of the full numerical analysis. 

Open with DEXTER  
In the text 
Fig. 8 γ(M) for NFW halos with masses ranging from 10^{9} M_{⊙}h^{1} to 10^{15} M_{⊙}h^{1} using Eq. (52) with . The four bands of lines show the results (from top to bottom) of ζ = 0.1,1.0,10,100, and within each band, we have (from top to bottom) n = 0.5,0.4,0.3, and 0.2. 

Open with DEXTER  
In the text 
Fig. 9 Rescaling parameter μ_{200} versus the maximum gravitational enhancement γ_{max} for several discussed models (solid colored lines) and the respective partially screened regions (colored semitransparent areas). The red line shows the space spanned by the HuSawicki f(R) model when varying  f_{R0}  from 10^{8} (left end) to 10^{4} (right end), while leaving n = 1 fixed. For the Symmetron, we also varied one parameter per curve (in blue, green, and violet), while leaving the others fixed at β = 1, L = 1 Mpc h^{1}, a_{ssb} = 1/2. The orange line indicates the region for nDGP models with r_{c} ∈ [0.5, 10] (left to right) with the associated semitransparent area spanning wider than the plot range. See Sect. 4.1 for details. 

Open with DEXTER  
In the text 