© 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 well-tested 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 screened-modified 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 self-similar 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 N-body 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.3, 0.7), ${\mathit{M}}_{\mathrm{Pl}}^{-2}\mathrm{\equiv }\mathrm{8}\mathit{\pi G}$, ${\mathit{\rho }}_{\mathrm{c}}\mathrm{=}\mathrm{3}{\mathit{H}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}$, $\overline{{\mathit{\rho }}_{\mathrm{m}}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{\rho }}_{\mathrm{c}}$ 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 Navarro-Frenk-White (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 Einstein-Hilbert action. Specifically, in scalar-tensor theories one (or multiple) additional scalar field is added with their respective Lagrangian L_ϕ, and the general action reads as $\begin{array}{ccc}\mathit{S}\mathrm{=}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{4}}\mathit{x}\sqrt{\mathrm{-}\mathit{g}}\left(\frac{{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}}{\mathrm{2}}\mathit{R}\mathrm{+}{\mathit{L}}_{\mathit{\varphi }}\mathrm{\left(}\mathit{\varphi ,\partial \varphi ,}{\mathit{\partial }}^{\mathrm{2}}\mathit{\varphi }\mathrm{\right)}\right)\mathrm{+}{\mathrm{𝒮}}_{\mathrm{m}}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{,}\begin{array}{c}\mathrm{˜}\\ {\mathit{g}}_{\mathit{\mu \nu }}\end{array}\mathrm{\right)}\mathit{,}& & \end{array}$(1)where the matter fields ψ⁽ⁱ⁾ couple to the metric $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{g}}_{\mathit{\mu \nu }}\end{array}\mathrm{\equiv }{\mathit{A}}^{\mathrm{2}}{\mathit{g}}_{\mathit{\mu \nu }}$, and R is the Ricci scalar. We limit ourselves to one scalar field and a conformal coupling to matter¹. Since the scalar field is coupled to matter, i.e., A = A(ϕ), an extra fifth force arises that is in the non-relativistic limit, and per unit mass is given by ${{F}}_{\mathit{\varphi }}\mathrm{=}\mathrm{-}{\nabla }\log \mathit{A}\mathit{.}$(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 ∂²ϕ. In this case, the screening is caused by the relation between the kinetic term of L_ϕ, i.e., $\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{\partial \varphi }{\mathrm{)}}^{\mathrm{2}}$, 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 ∂²ϕ, and the k-mouflage mechanism (Babichev et al. 2009; Brax & Valageas 2014a) , where the magnitude of (∂ϕ)² 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 N-body simulations. For theories that screen via the first method described above (chameleon-like 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 environment-dependent dilation model (Brax et al. 2011, 2012c). For models that employ the second screening method (Vainshtein-like 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 k-moflage (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 $\mathit{\gamma }\mathrm{\equiv }{\mathit{G}}_{\mathrm{eff}}\mathit{/}\mathit{G}\mathrm{-}\mathrm{1}$(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.

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 $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{{\mathit{\varphi }}^{\mathrm{2}}}{{\mathit{M}}^{\mathrm{2}}}\\ \mathit{V}\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\mu }}^{\mathrm{2}}{\mathit{\varphi }}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}\mathit{\lambda }{\mathit{\varphi }}^{\mathrm{4}}\mathit{.}\end{array}$Here, the prefactors M, μ and λ can be rewritten as $\begin{array}{ccc}& & {\mathit{M}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}{\mathrm{\Omega }}_{\mathrm{m}\mathrm{0}}{\mathit{\rho }}_{\mathrm{c}\mathrm{0}}{\mathit{L}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}}\mathit{,}\\ & & \mathit{\mu }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}}\mathit{L}}\mathit{,}\\ & & {\mathit{\lambda }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}{\mathit{M}}_{\mathrm{Pl}}}{\sqrt{\mathrm{8}}{\mathit{L}}^{\mathrm{3}}{\mathrm{\Omega }}_{\mathrm{m}\mathrm{0}}{\mathit{\rho }}_{\mathrm{c}\mathrm{0}}\mathit{\beta }}\mathit{,}\end{array}$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 Friedmann-Lemaître-Robertson-Walter (FLRW) background and in the quasi-static limit (Noller et al. 2014) , the field equation becomes ${\mathrm{\nabla }}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{\varphi }\end{array}\mathrm{=}\frac{{\mathit{a}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{2}}}\left(\begin{array}{c}\mathrm{˜}\\ \mathit{\varphi }\end{array}\frac{{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}}{{\mathit{a}}^{\mathrm{3}}}\frac{{\mathit{\rho }}_{\mathrm{m}}}{\overline{{\mathit{\rho }}_{\mathrm{m}}}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{\varphi }\end{array}\mathrm{+}\stackrel{}{{{\mathit{\varphi }}^{\mathrm{3}}}_{\mathrm{˜}}}\right)$(9)where $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\varphi }\end{array}\mathrm{=}\sqrt{\mathit{\lambda }}\mathit{\varphi }\mathit{/}\mathit{\mu }$. The fifth force is given by $\begin{array}{ccc}{\mathit{F}}_{\mathit{\varphi }}\mathrm{=}\frac{\mathit{\varphi }\mathrm{\nabla }\mathit{\varphi }}{{\mathit{M}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(10)In regions of high ambient matter density (${\mathit{\rho }}_{\mathrm{m}}\mathrm{\gg }{\mathit{\mu }}^{\mathrm{2}}{\mathit{M}}^{\mathrm{2}}\mathrm{=}{\mathit{\rho }}_{\mathrm{c}\mathrm{0}}\mathit{/}{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}$), the field will reside close to φ = 0, and the fifth force will be screened. In low-density 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 ${\mathit{\gamma }}_{\max }\mathrm{=}\mathrm{2}{\mathit{\beta }}^{\mathrm{2}}\left[\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{a}}_{\mathrm{ssb}}}{\mathit{a}}\right)}^{\mathrm{3}}\right]\mathit{.}$(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. Hu-Sawicki f(R) model

The f(R) gravity models are a group of modified gravity theories where the Ricci scalar R in the Einstein-Hilbert 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² translates f(R)-gravity into a scalar tensor theory. The theory screens via the chameleon mechanism.

The fifth force is given by ${\mathit{F}}_{\mathit{\varphi }}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\nabla }{\mathit{f}}_{\mathit{R}}$(12)where f_R ≡ df(R)/ dR is the scalar field. In a FLRW background spacetime and in the quasi-static limit (Noller et al. 2014) , the field-equation becomes $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{f}}_{\mathit{R}}& \mathrm{=}& \mathrm{-}{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\left(\frac{{\mathit{\rho }}_{\mathrm{m}}}{\overline{{\mathit{\rho }}_{\mathrm{m}}}}\mathrm{-}\mathrm{1}\right)\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}\\ & & \mathrm{\times }\left[\left(\mathrm{1}\mathrm{+}\mathrm{4}\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_{\mathrm{m}}}\right){\left(\frac{{\mathit{f}}_{\mathit{R}\mathrm{0}}}{{\mathit{f}}_{\mathit{R}}}\right)}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{a}}^{-3}\mathrm{-}\mathrm{4}\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_{\mathrm{m}}}\right]\mathit{,}\end{array}$(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 ${\mathit{\gamma }}_{\max }\mathrm{=}\mathrm{2}{\mathit{\beta }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathit{,}$(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 four-dimensional brane that itself is embedded in a five-dimensional spacetime. The DGP model can give rise to self-acceleration 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 $\begin{array}{ccc}{\mathit{F}}_{\mathit{\varphi }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\nabla }\mathit{\varphi ,}& & \end{array}$(15)where the evolution of the scalar field ϕ, the co-called brane-bending mode, is determined by the field equation³$\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\varphi }\mathrm{+}\frac{{\mathit{r}}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{a}}^{\mathrm{2}}}\mathrm{(}\mathrm{(}\mathrm{\nabla }\mathit{\varphi }{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathrm{\nabla }}_{\mathit{i}}{\mathrm{\nabla }}_{\mathit{j}}\mathit{\varphi }{\mathrm{)}}^{\mathrm{2}}\mathrm{)}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{a}{\mathit{\beta }}_{\mathrm{DGP}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{\delta }}_{\mathrm{m}}\mathit{,}& & \end{array}$(16)where $\begin{array}{ccc}{\mathit{\beta }}_{\mathrm{DGP}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{2}\mathrm{\left(}{\mathit{r}}_{\mathrm{c}}{\mathit{H}}_{\mathrm{0}}\mathrm{\right)}\frac{\mathit{H}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{\mathit{H}}_{\mathrm{0}}}\left(\mathrm{1}\mathrm{+}\frac{\mathit{Ḣ}}{\mathrm{3}{\mathit{H}}^{\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(17)The maximum value of the gravitational force in this model is $\begin{array}{ccc}{\mathit{\gamma }}_{\max }\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}}\mathrm{·}& & \end{array}$(18)For high values of r_c, we have β_DGP(a = 1) ∝ r_c, so ${\mathit{\gamma }}_{\max }\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{c}}}\mathit{.}$ 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 $\begin{array}{ccc}\frac{{\mathit{\varphi }}^{\mathrm{\prime }}}{\mathit{r}}\mathrm{+}\frac{\mathrm{2}{\mathit{r}}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{a}}^{\mathrm{2}}}{\left(\frac{{\mathit{\varphi }}^{\mathrm{\prime }}}{\mathit{r}}\right)}^{\mathrm{2}}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{a}{\mathit{\beta }}_{\mathrm{DGP}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathit{\delta }}_{\mathrm{m}}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{r}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{·}& & \end{array}$(19)This translates to a solution for the fifth force profile given by $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}\frac{\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varphi }}^{\mathrm{\prime }}}{{\mathrm{\Phi }}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}}\left[\frac{\mathrm{2}\mathrm{(}\mathrm{-}\mathrm{1}\mathrm{+}\sqrt{\mathrm{1}\mathrm{+}\mathit{ϵ}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{)}}{\mathit{ϵ}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right]& & \end{array}$(20)where $\begin{array}{ccc}\mathit{ϵ}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{8}\mathrm{(}{\mathit{r}}_{\mathrm{c}}{\mathit{H}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{a}}^{\mathrm{3}}}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathit{\delta }}_{\mathrm{m}}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{r}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{·}& & \end{array}$(21)We have two regimes. For ϵ ≪ 1 (large r), we get $\mathit{\gamma }\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}}$ and gravity is maximally modified. On the other hand, for ϵ ≫ 1 (small r), we have $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}}\frac{\mathrm{2}}{\sqrt{\mathit{ϵ}}}\mathrm{\ll }\frac{\mathrm{1}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}}\mathit{,}& & \end{array}$(22)and the fifth-force is screened.

2.3. The unified { m(a),β(a) } description of chameleon-like models

Brax et al. (2012b,a) show that any scalar-tensor 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) $\begin{array}{ccc}& & \mathit{V}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathit{V}}_{\mathrm{0}}\mathrm{-}\frac{\mathrm{3}}{{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}}{\mathrm{\int }}_{{\mathit{a}}_{\mathit{i}}}^{\mathit{a}}\frac{{\mathit{\beta }}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}{\mathit{m}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{\rho }}_{\mathrm{m}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{a}\\ & & \mathit{\varphi }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathit{\varphi }}_{\mathit{i}}\mathrm{+}\frac{\mathrm{3}}{{\mathit{M}}_{\mathrm{Pl}}}{\mathrm{\int }}_{{\mathit{a}}_{\mathit{i}}}^{\mathit{a}}\frac{\mathit{\beta }\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}{\mathit{m}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{\rho }}_{\mathrm{m}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{a}\\ & & \log \mathit{A}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\log \mathit{A}\mathrm{\left(}{\mathit{a}}_{\mathit{i}}\mathrm{\right)}\mathrm{+}\frac{\mathrm{3}}{{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}}{\mathrm{\int }}_{{\mathit{a}}_{\mathit{i}}}^{\mathit{a}}\frac{{\mathit{\beta }}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}{\mathit{m}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{\rho }}_{\mathrm{m}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a,}\end{array}$where ${\mathit{\rho }}_{\mathrm{m}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{3}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}\mathit{/}{\mathit{a}}^{\mathrm{3}}$.

One of the advantages of using this form is the direct relationship between the evolution of linear matter perturbations and the β, m functions $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ {\mathit{\delta }}_{\mathrm{m}}\end{array}\mathrm{+}\mathrm{2}\mathit{H}\mathit{\delta ̇}\mathrm{m}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{a}}^{-3}{\mathit{\delta }}_{\mathrm{m}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{2}{\mathit{\beta }}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{k}}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\mathit{m}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}\right)\mathrm{·}& & \end{array}$(26)The two functions are, in the case of the Symmetron model, given by $\begin{array}{ccc}\frac{\mathit{m}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{=}& \frac{\mathit{\beta }\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{\mathit{\beta }}_{\mathrm{0}}}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathrm{(}{\frac{{\mathit{a}}_{\mathrm{ssb}}}{\mathit{a}}}^{\mathrm{)}}}^{\mathrm{3}}}\mathrm{;}& \end{array}$(27)i.e., they follow the same evolution. The parameters β₀ and m₀ are related to the model parameters described in the previous section via β₀ = β and ${\mathit{m}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathit{L}}$.

For the Hu-Sawicki f(R) model, on the other hand, we have $\begin{array}{ccc}\mathit{\beta }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{6}}}& & \end{array}$(28)and $\begin{array}{ccc}\mathit{m}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\frac{{\mathit{H}}_{\mathrm{0}}\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{+}\mathrm{4}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}{\sqrt{\mathrm{\left|}{\mathit{f}}_{\mathit{R}\mathrm{0}}\mathrm{\right|}}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{a}}^{-3}\mathrm{+}\mathrm{4}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{+}\mathrm{4}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}\right)}^{\mathit{n}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{2}}\mathit{,}& & \end{array}$(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/$\sqrt{\mathrm{6}{\mathit{\beta }}_{\mathrm{DGP}}}$.

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 ${\mathit{a}}_{\mathit{\rho }}\mathrm{=}{\mathit{a}}_{\mathrm{env}}\mathrm{(}\overline{{\mathit{\rho }}_{\mathrm{m}}}\mathit{/}\mathit{\rho }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$. The screening condition for an object with Newtonian potential Φ_N in a region of density ρ can therefore be written as⁴$\begin{array}{ccc}\mathit{ϵ}\mathrm{\left(}{\mathit{a}}_{\mathit{\rho }}\mathit{,}{\mathit{a}}_{\mathrm{env}}\mathrm{\right)}\mathrm{=}\min \left[\frac{\mathrm{\left|}\mathit{\varphi }\mathrm{\right(}{\mathit{a}}_{\mathit{\rho }}\mathrm{)}\mathrm{-}\mathit{\varphi }\mathrm{(}{\mathit{a}}_{\mathrm{env}}\mathrm{\left)}\mathrm{\right|}}{\mathrm{2}\mathit{\beta }\mathrm{\left(}{\mathit{a}}_{\mathrm{env}}\mathrm{\right)}{\mathit{M}}_{\mathrm{Pl}}{\mathrm{\Phi }}_{\mathrm{N}}}\mathit{,}\mathrm{1}\right]\mathit{.}& & \end{array}$(30)An object is screened whenever ϵ(a_ρ,a_env) ≪ 1. Here, ${\mathit{a}}_{\mathrm{env}}\mathrm{=}\mathit{a}\mathrm{(}\overline{{\mathit{\rho }}_{\mathrm{m}}}\mathit{/}{\mathit{\rho }}_{\mathrm{env}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$ 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 $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}\mathrm{2}{\mathit{\beta }}^{\mathrm{2}}\mathrm{\left(}{\mathit{a}}_{\mathrm{env}}\mathrm{\right)}\mathit{ϵ}\mathrm{\left(}{\mathit{a}}_{\mathit{\rho }}\mathit{,}{\mathit{a}}_{\mathrm{env}}\mathrm{\right)}\mathit{.}& & \end{array}$(31)We should note that this is a very rough analytical approximation, but it is usually good enough for order-of-magnitude estimates. To get accurate predictions, we need to solve the field equation numerically, as discussed in the next section.

Fig. 1 γ vs. R/R200 for NFW halos with masses ranging from 109 M⊙h-1 to 1014.5 M⊙h-1. The left panel shows the Symmetron solution using (assb,L,β) = (0.7, 1 Mpc h-1, 1), the central panel the Hu & Sawickif(R) solution with (| fR0 | ,n) = (10-6, 1), and the right panel the nDGP solution with rc = 1.

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 NFW-equations), and then use Newton-Gauss Seidel relaxation with multigrid acceleration to obtain the solution.

In spherical symmetry the field equation becomes one-dimensional, φ = φ(r) and ${\mathrm{\nabla }}^{\mathrm{2}}\mathit{\phi }\mathrm{=}\frac{{\mathit{d}}^{\mathrm{2}}\mathit{\phi }}{\mathrm{d}{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}}{\mathit{r}}\frac{\mathrm{d}\mathit{\phi }}{\mathrm{d}\mathit{r}}$. For the Symmetron model, we implemented a simple fixed-grid multigrid solver for this one-dimensional problem. To do this we write Eq. (9) as ℒ_i = 0 where $\begin{array}{ccc}{\mathrm{ℒ}}_{\mathit{i}}\mathrm{=}\mathrm{(}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\varphi }{\mathrm{)}}_{\mathit{i}}\mathrm{-}\frac{{\mathit{a}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{2}}}\left({\mathit{\varphi }}_{\mathit{i}}\frac{{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}}{{\mathit{a}}^{\mathrm{3}}}\frac{{\mathit{\rho }}_{\mathrm{m}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}}{\overline{{\mathit{\rho }}_{\mathrm{m}}}}\mathrm{-}{\mathit{\varphi }}_{\mathit{i}}\mathrm{+}{\mathit{\varphi }}_{\mathit{i}}^{\mathrm{3}}\right)\mathit{,}& & \end{array}$(32)where subscript i denotes the value at the ith gridpoint and $\begin{array}{ccc}\mathrm{(}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\varphi }{\mathrm{)}}_{\mathit{i}}\mathrm{=}\frac{{\mathit{\varphi }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{+}{\mathit{\varphi }}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{-}\mathrm{2}{\mathit{\varphi }}_{\mathit{i}}}{\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}}{{\mathit{r}}_{\mathit{i}}}\frac{{\mathit{\varphi }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{\varphi }}_{\mathit{i}\mathrm{-}\mathrm{1}}}{\mathrm{\Delta }\mathit{r}}& & \end{array}$(33)is a second-order discretization of the Laplacian. We start by making a guess for the solution, and then we loop through the grid updating the solution using $\begin{array}{ccc}{\mathit{\varphi }}_{\mathit{i}}\mathrm{\leftarrow }{\mathit{\varphi }}_{\mathit{i}}\mathrm{-}\frac{{\mathrm{ℒ}}_{\mathit{i}}}{\mathit{\partial }{\mathrm{ℒ}}_{\mathit{i}}\mathit{/}\mathit{\partial }{\mathit{\varphi }}_{\mathit{i}}}\mathrm{·}& & \end{array}$(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 Semi-analytical solutions for γ for the Symmetron (left) and f(R) (right) using the same parameter values (and color coding) as in Fig. 1.

The f(R) equation Eq. (13) is stiffer than the Symmetron equation for low values of | f_R0 |, and a simple fixed-grid one-dimensional 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 well-tested 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³ 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⁵. 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 semi-analytically in the { m(a),β(a) } formulations in Sect. 3.4.

Fig. 3 vs halo mass M. The left panel shows the Symmetron curves with (assb,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 | fR0 | = (−4, −5, −6, −7) (in red, blue, green, and purple, respectively) and the right panel the nDGP solution with rc = (0.5, 1 ,3) (from top to bottom).

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₂₀₀ and drops again for larger radii. This drop happens due to the fifth force having a finite range, ${\mathit{F}}_{\mathit{\phi }}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\mathit{mr}}$.

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¹² M_⊙ seem to be fully screened in the central (R ≲ 0.1R₂₀₀) 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 semi-analytical 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₂₀₀ where γ peaks. The semi-analytic 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 concentration-mass 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⁶.

3.2. Profile variation with halo mass

To capture the change of the fifth force profile with halo mass, we chose to show the mass-weighted mean of the fifth force $\mathit{\gamma ̅}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{x}}_{\mathrm{cut}}}\hspace{0.17em}\mathrm{d}\mathit{x\gamma }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{w}\mathrm{\left(}\mathit{x}\mathrm{\right)}$(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 one-tenth 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 one-third. For lighter halos falls off rapidly. In the righter-most 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 mass-concentration 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 μ₂₀₀ in the following way: $\mathit{\gamma ̅}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{200}}\mathrm{\right)}\mathrm{\equiv }\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\gamma }}_{\max }$(36)where is still the mass-weighted average of the enhancement of the fifth force as defined in Eq. (35) In general μ₂₀₀ is a function of the various model parameters. Halos with mass much higher than μ₂₀₀ are in the fully screened regime, whereas halos with mass much lower than μ₂₀₀ are in the unscreened regime.

In this section we find the functional form of μ₂₀₀ for each model merely empirically from our numerical solutions. However, in the next section, we derive semi-analytic 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 μ₂₀₀. For example, for the f(R) model with (n, | f_R0 | ) = (1, 10^-6) (in green), γ( ~10¹³ M_⊙ h^-1) = 1/6 = 1/2γ_max, i.e., μ₂₀₀, is in that particular case ~10¹³ M_⊙ h^-1.

For the Symmetron model, we found μ₂₀₀ to be well fit by ${\mathit{\mu }}_{\mathrm{200}}^{\mathrm{Symmetron}}\mathrm{=}\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{10}}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{h}}^{-1}\mathrm{\times }{\left(\frac{\mathit{L}}{\hspace{0.17em}\mathrm{Mpc}\hspace{0.17em}{\mathit{h}}^{-1}}\right)}^{\mathrm{3}}{\mathit{a}}_{\mathrm{ssb}}^{-4.5}\mathit{,}$(37)and γ_max is given by Eq. (11). Figure 4 shows the resulting rescaled curves, i.e., versus M/μ₂₀₀. 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 ≪ μ₂₀₀ through γ(M ≈ μ₂₀₀) = γ_max/ 2 to γ ~ 0 for M ≫ μ₂₀₀ – 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.

Fig. 5 versus M/μ200 for the Hu-Sawicki 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 fR0. μ200 for this model is given by Eq. (38).

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 ${\mathit{\mu }}_{\mathrm{200}}^{\mathit{f}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{13}}{\mathit{M}}_{\mathrm{\odot }}{\mathrm{h}}^{-1}\mathrm{\times }{\left(\frac{\mathrm{\left|}{\mathit{f}}_{\mathit{R}\mathrm{0}}\mathrm{\right|}}{{\mathrm{10}}^{-6}}\right)}^{\mathrm{1.5}}$(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 color-coded 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 ≪ μ₂₀₀, γ ≈ γ_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 μ₂₀₀ 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 ${\log }_{\mathrm{10}}{\mathit{\mu }}_{\mathrm{200}}\mathrm{=}\mathit{A}\mathrm{-}\mathrm{(}{\mathit{r}}_{\mathrm{c}}\mathit{/}\mathit{a}{\mathrm{)}}^{\mathit{b}}{\mathrm{e}}^{\mathrm{-}{\mathit{r}}_{\mathrm{c}}\mathit{/}\mathit{a}}\mathit{,}$(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.

Figure 6 once more shows the difference between the Symmetron and f(R) model with their characteristic “inverse-S-shape” on the one side, and, the nDGP – with a log-linear relation – on the other. A simple fitting function that is able to match the characteristic “inverse-S-shape” given by $\begin{array}{ccc}\mathit{\gamma ̅}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{=}{\mathit{\gamma }}_{\max }\left[\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{\mathit{\pi }}\arctan \mathrm{(}{\mathit{M}}_{\mathrm{200}}\mathit{/}{\mathit{\mu }}_{\mathrm{200}}\mathrm{)}\right]\mathit{.}& & \end{array}$(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 $\begin{array}{ccc}\mathit{\gamma ̅}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{=}\mathrm{𝒜}{\mathit{\gamma }}_{\max }\left(\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{\mathit{\pi }}\arctan \left[\mathrm{(}{\mathit{M}}_{\mathrm{200}}\mathit{/}{\mathit{\mu }}_{\mathrm{200}}{\mathrm{)}}^{\mathit{\alpha }}\right]\right)& & \end{array}$(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 “inverse-S-shape” flips vertically and turns into an “S-shape”. Quantitatively, we can relate α to the half width of the partially screened region through $\mathit{W}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\tan {\left(\frac{\mathit{\pi d}}{\mathrm{2}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{\left|}\mathit{\alpha }\mathrm{\right|}}\mathrm{·}$(42)Here, M₂₀₀ ∈ [μ₂₀₀/W,Wμ₂₀₀] 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, two-parameter 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. Semi-analytical derivation of μ₂₀₀ 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 μ₂₀₀. It turns out that we can do this quite generally, so our results below can be applied to any scalar-tensor theory with a Lagrangian $\mathit{L}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{\partial \varphi }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathit{V}\mathrm{(}\mathit{\varphi }\mathrm{)}$ 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 semi-analytic approximation, which formally holds only for a spherical top-hat overdensity, is valid for NFW halos. In Fig. 2 we show the semi-analytical 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 mass-weighted average of γ, which reads as $\begin{array}{ccc}\mathit{\gamma ̅}\mathrm{\simeq }{\mathit{\gamma }}_{\max }\frac{\mathit{M}\mathrm{(}\mathit{<}{\mathit{x}}_{\mathrm{screen}}\mathrm{)}}{\mathit{M}\mathrm{(}\mathit{<}{\mathit{x}}_{\mathrm{cut}}\mathrm{)}}\mathrm{·}& & \end{array}$(43)Furthermore, using the definition of μ₂₀₀, namely $\mathit{\gamma ̅}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{200}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\gamma }}_{\max }$, we require $\begin{array}{ccc}\frac{\mathit{M}\mathrm{(}\mathit{<}{\mathit{x}}_{\mathrm{screen}}\mathrm{)}}{\mathit{M}\mathrm{(}\mathit{<}{\mathit{x}}_{\mathrm{cut}}\mathrm{)}}{|}_{{\mathit{M}}_{\mathrm{200}}\mathrm{=}{\mathit{\mu }}_{\mathrm{200}}}\mathrm{\simeq }\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}& & \end{array}$(44)which requires for typical halo masses.

The dependence of x_screen, hence of μ₂₀₀, on the various SMG model parameters can be calculated using the screening condition Eq. (30). In this case, the condition reads as $\begin{array}{ccc}\frac{\mathrm{\left|}\mathit{\varphi }\mathrm{\right(}{\mathit{a}}_{\mathit{\rho }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{)}\mathrm{-}\mathit{\varphi }\mathrm{(}{\mathit{a}}_{\mathrm{env}}\mathrm{\left)}\mathrm{\right|}}{\mathrm{2}\mathit{\beta }\mathrm{\left(}{\mathit{a}}_{\mathrm{env}}\mathrm{\right)}{\mathit{M}}_{\mathrm{Pl}}{\mathrm{\Phi }}_{\mathrm{N}}\mathrm{\left(}{\mathit{x}}_{\mathrm{screen}}\mathrm{\right)}}\mathrm{=}\mathrm{1}& & \\ {\mathit{a}}_{\mathit{\rho }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{a}}_{\mathrm{env}}{\left(\frac{\overline{{\mathit{\rho }}_{\mathrm{m}}}}{\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathrm{screen}}\mathrm{\right)}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}& & \end{array}$With ${\mathrm{\Phi }}_{\mathrm{N}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{GM}\mathrm{(}\mathit{<}\mathit{r}\mathrm{)}}{\mathit{r}}$ and M₂₀₀ = 4π/ 3Δ_virΩ_mρ_c0r³ → r ∝ M^1/3, we can easily calculate how μ₂₀₀ scales with the model parameters in any given model.

For the Symmetron we find the scaling M( <${\mathit{x}}_{\mathrm{screen}}{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{\propto }\frac{{\mathit{L}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}}$. As a result, μ₂₀₀ scales as $\begin{array}{ccc}{\mathit{\mu }}_{\mathrm{200}}^{\mathrm{Symmetron}}\mathrm{=}\mathrm{(}\mathrm{~}\mathrm{0.4}\mathrm{\pm }\mathrm{0.3}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{10}}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{h}}^{-1}\mathrm{\times }{\left(\frac{\mathit{L}}{{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{Mpc}\mathit{/}\mathit{h}}\right)}^{\mathrm{3}}& & \end{array}$(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 Semi-analytical predictions for versus M/μ200 for the Hu-Sawicki 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.

For the Hu-Sawicki f(R) model, we find the scaling μ₂₀₀ ∝ | f_R0 | ^3/2C^{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 μ₂₀₀. Again the proportionality factor is found numerically and gives us $\begin{array}{ccc}{\mathit{\mu }}_{\mathrm{200}}^{\mathit{f}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\mathrm{=}\mathrm{(}\mathrm{~}\mathrm{1.1}\mathrm{\pm }\mathrm{0.4}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{13}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}{\mathit{h}}^{-1}\mathrm{\times }{\left(\frac{\mathrm{\left|}{\mathit{f}}_{\mathit{R}\mathrm{0}}\mathrm{\right|}}{{\mathrm{10}}^{-6}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(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 $\mathit{\phi }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{~}\frac{\mathit{\beta }\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{\mathit{a}}^{\mathrm{3}}{\mathit{m}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}}$. This in the condition Eq. (45) means that we generally expect ${\mathit{\mu }}_{\mathrm{200}}\mathrm{\propto }{\mathit{\lambda }}_{\mathit{\phi }}^{\mathrm{3}}\mathrm{\left(}\mathit{a}\mathrm{\right)}$ where λ_φ(a) = m^-1(a) is the range of the fifth force.

In Fig. 7 we show the semi-analytical profiles for in terms of M₂₀₀/μ₂₀₀ 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 $\begin{array}{ccc}\mathit{m}\mathrm{\left(}\mathit{a}\mathrm{\right)}& \mathrm{=}& {\mathit{m}}_{\mathrm{0}}{\mathit{a}}^{\mathrm{-}\mathit{r}}\\ \mathit{\beta }\mathrm{\left(}\mathit{a}\mathrm{\right)}& \mathrm{=}& {\mathit{\beta }}_{\mathrm{0}}{\mathrm{e}}^{\frac{\mathit{s}}{\mathrm{2}\mathit{r}\mathrm{-}\mathrm{3}}\mathrm{(}{\mathit{a}}^{\mathrm{2}\mathit{r}\mathrm{-}\mathrm{3}}\mathrm{-}\mathrm{1}\mathrm{)}}\end{array}$where r,s,m₀, and β₀ 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 semi-analytically, so we leave this for future work.

Fig. 8 γ(M) for NFW halos with masses ranging from 109 M⊙h-1 to 1015 M⊙h-1 using Eq. (52) with $\mathit{G}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{\left(}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{\zeta z}{\mathrm{)}}^{\mathit{n}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{\zeta z}}$. 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.

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), second-order equations of motion the field equation in spherical symmetry are integrable, leading to an algebraic equation on the form $\begin{array}{ccc}\mathit{f}\mathrm{(}{\mathit{\phi }}^{\mathrm{\prime }}\mathit{/}\mathit{r}\mathrm{)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathit{\delta }}_{\mathrm{m}}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{=}\frac{{\mathrm{\Delta }}_{\mathrm{vir}}}{\mathrm{3}{\mathit{x}}^{\mathrm{3}}}\left(\frac{\log \mathrm{(}\mathrm{1}\mathrm{+}\mathit{cx}\mathrm{)}\mathrm{-}\frac{\mathit{cx}}{\mathrm{1}\mathrm{+}\mathit{cx}}}{\log \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\frac{\mathit{c}}{\mathrm{1}\mathrm{+}\mathit{c}}}\right)& & \end{array}$(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 $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}{\mathit{\gamma }}_{\max }\mathit{G}\left(\frac{{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}\mathrm{\right)}}{{\mathit{x}}^{\mathrm{3}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{cx}\mathrm{\right)}}\right)& & \end{array}$(52)for some function G that satisfies G(0) = 1 and G(∞) = 0. For DGP we have $\mathit{G}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{(}\sqrt{\mathrm{1}\mathrm{+}\mathit{\zeta z}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{\zeta z}}$ with $\mathit{\zeta }\mathrm{=}\frac{\mathrm{8}{\mathit{r}}_{\mathrm{c}}^{\mathrm{2}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}{\mathrm{\Delta }}_{\mathrm{vir}}}{\mathrm{9}{\mathit{\beta }}_{\mathrm{DGP}}^{\mathrm{2}}}$.

For models that only have a shift symmetry, like the k-mouflage screening mechanism, the field equation, $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left[{\mathit{r}}^{\mathrm{2}}\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}\frac{\mathrm{d}\mathit{\phi }}{\mathrm{d}\mathit{r}}\right]\mathrm{=}\frac{\mathit{\beta \delta }{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{M}}_{\mathrm{Pl}}}\mathit{,}& & \end{array}$(53)can also be integrated up, and from this it follows that the fifth-force F_φ ∝ ∇φ is a function of $\frac{{{\mathrm{\int }}^{}}_{\mathrm{0}}^{\mathit{r}}{\mathit{\delta }}_{\mathrm{m}}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}}{{\mathit{r}}^{\mathrm{2}}}$ for spherical symmetry. This leads to a slightly different form $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}{\mathit{\gamma }}_{\max }\mathit{F}\left(\frac{{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}\mathrm{\right)}}{{\mathit{x}}^{\mathrm{2}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{cx}\mathrm{\right)}}\right)& & \end{array}$(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 $\mathit{G}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\zeta z}{\mathrm{)}}^{\mathit{n}}\mathrm{-}\mathrm{1}}{\mathit{\zeta z}}$ 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 semi-transparent areas). The red line shows the space spanned by the Hu-Sawicki f(R) model when varying | fR0 | 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, assb = 1/2. The orange line indicates the region for nDGP models with rc ∈ [0.5, 10] (left to right) with the associated semi-transparent area spanning wider than the plot range. See Sect. 4.1 for details.

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 N-body 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 μ₂₀₀ (Sect. 3.3) demonstrates that in spite of their foundational differences, the three screened-modified gravity models investigated show self-similar 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 μ₂₀₀ 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,μ₂₀₀) ≈ (1.75, 4.5 × 10¹¹ M_⊙h^-1). The semi-transparent 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 μ₂₀₀, the semi-transparent 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, μ₂₀₀/W can be interpreted as the minimum mass where one expects a variation in halo properties; that is to say:

• For M₂₀₀ ≲ μ₂₀₀/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₂₀₀ ∈ [μ₂₀₀/W,Wμ₂₀₀] 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 N-body simulations). Also, comparing the inner and outer regions of individual objects can lead to insight (for instance, the ratio M₅₀₀/M₂₀₀). 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μ₂₀₀, the halo is fully screened, so no deviation from the ΛCDM predictions is expected.

Conclusively, in any observable-mass scaling relation, we predict a peculiar feature at μ₂₀₀, 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 screened-modified gravity considered, but differs significantly between them. In particular, the value of α found in the numerical fits of Sect. 3.3 translates to log ₁₀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 as-sumed spherical symmetry and perfect NFW halos.However, in reality clusters of galaxies do not exactlyfollow an NFW mass distribution, and halos also pos-sess a nonzero ellipticity (e.g., Oguriet al. 2010). Another phenomenonfalling in the same category is our assignment of a unique mass-concentration 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, individ-ual clusters (as any astrophysical object) have mostly a limitedexplanatory power and must be stacked for analysis.

• Model assumptions. To rescale the curves accordingly, an intermediate-mass 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 μ₂₀₀ 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 N-body 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) $\begin{array}{ccc}\mathit{\gamma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{Dyn}}\mathrm{(}\mathit{<}\mathit{r}\mathrm{)}}{{\mathit{M}}_{\mathrm{Lensing}}\mathrm{(}\mathit{<}\mathit{r}\mathrm{)}}\mathrm{-}\mathrm{1}\mathit{,}& & \end{array}$(55)where $\begin{array}{ccc}& & {\mathit{M}}_{\mathrm{Dyn}}\mathrm{\equiv }\mathrm{4}\mathit{\pi }\mathrm{\int }{\mathit{r}}^{\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}\left(\frac{\mathrm{\Phi }\mathrm{+}\mathrm{\Psi }}{\mathrm{2}}\right)\hspace{0.17em}\mathrm{d}\mathit{r}\\ & & {\mathit{M}}_{\mathrm{Lensing}}\mathrm{\equiv }\mathrm{4}\mathit{\pi }\mathrm{\int }{\mathit{r}}^{\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}\left(\frac{\mathrm{\Phi }\mathrm{-}\mathrm{\Psi }}{\mathrm{2}}\right)\hspace{0.17em}\mathrm{d}\mathit{r}\mathit{.}\end{array}$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⁷. 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 power-law relation obtained from (ΛCDM) N-body simulations (e.g., Becker et al. 2007). The lensing mass can be read from the stacked density profiles assuming a constant mass-richness relation (Johnston et al. 2007; Hoekstra et al. 2013). This means that our mass-weighted 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)⁸.

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 Hu-Sawicki f(R), and the nDGP models (Sect. 3.3). In addition, we found a semi-analytic 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 Hu-Sawicki f(R), and the nDGP models with three parameters: (i) the center of the partially screened regime μ₂₀₀; (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 Hu-Sawicki f(R), and the Dilaton models.

Quantifying observational constraints in the proposed (μ₂₀₀,γ_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.

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

Appendix A: NFW equations

The Navarro-Frenk-White (NFW) profile is an empirically derived halo profile (Navarro et al. 1995, 1996; 1997) from N-body simulations. The NFW density profile reads as a function of normalized distance from the center x ≡ r/R_vir: $\frac{\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}}{{\mathit{\rho }}_{\mathrm{c}\mathrm{0}}}\mathrm{=}\frac{{\mathrm{\Delta }}_{\mathrm{vir}}{\mathit{c}}^{\mathrm{2}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}\mathrm{\right)}}{\mathrm{3}\mathit{x}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{cx}{\mathrm{)}}^{\mathrm{2}}}$(A.1)where ${\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\mathit{c}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}}\mathrm{·}$(A.2)Here, c is the concentration parameter, which we link to M₂₀₀ through the observed relationship found in Bullock et al. (2001): $\mathit{c}\mathrm{=}\mathrm{9}\mathrm{\times }{\left(\frac{{\mathit{M}}_{\mathrm{200}}}{\mathrm{3.2}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathit{h}}^{-1}}\right)}^{-0.13}\mathrm{·}$(A.3)Throughout this work we define Δ_vir ≡ 200, which means R_vir = R₂₀₀ and ${\mathit{M}}_{\mathrm{vir}}\mathrm{=}{\mathit{M}}_{\mathrm{200}}\mathrm{=}\mathrm{200}\mathrm{\times }\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi }{\mathit{R}}_{\mathrm{vir}}^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{c}\mathrm{0}}$.

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 $\frac{\hspace{0.17em}\mathrm{d}\mathit{M}}{\mathit{M}\mathrm{(}\mathit{<}{\mathit{x}}_{\mathrm{cut}}\mathrm{)}}\mathrm{=}\frac{{\mathit{c}}^{\mathrm{2}}\mathit{x}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}{\mathit{x}}_{\mathrm{cut}}\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{cx}{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{x,}$(A.4)where x_cut is defined below Eq. (35). This will be useful when considering mass-weighted quantities.

To calculate the screening condition we also need the expression for Φ_N. For a NFW halo, we find $\mathrm{\left|}{\mathrm{\Phi }}_{\mathrm{N}}\mathrm{\right|}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{\mathit{a}}\frac{\mathrm{\left(}{\mathit{H}}_{\mathrm{0}}{\mathit{R}}_{\mathrm{vir}}{\mathrm{)}}^{\mathrm{2}}{\mathrm{\Delta }}_{\mathrm{vir}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\right(}\mathit{c}\mathrm{)}}{\mathit{x}}\log \mathrm{(}\mathrm{1}\mathrm{+}\mathit{cx}\mathrm{)}\mathit{.}$(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 $\begin{array}{ccc}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathit{\delta }}_{\mathrm{m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{=}\frac{{\mathrm{\Delta }}_{\mathrm{vir}}}{\mathrm{3}{\mathit{x}}^{\mathrm{3}}}\left(\frac{\log \mathrm{(}\mathrm{1}\mathrm{+}\mathit{cx}\mathrm{)}\mathrm{-}\frac{\mathit{cx}}{\mathrm{1}\mathrm{+}\mathit{cx}}}{\log \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\frac{\mathit{c}}{\mathrm{1}\mathrm{+}\mathit{c}}}\right)\mathrm{=}\frac{{\mathrm{\Delta }}_{\mathrm{vir}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}\mathrm{\right)}}{{\mathit{x}}^{\mathrm{3}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{cx}\mathrm{\right)}}& & \end{array}$(B.1)where $\mathit{x}\mathrm{=}\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{vir}}}$. This leads to the fifth force profile being given by $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}{\mathit{\gamma }}_{\max }\left[\frac{\mathrm{2}\mathrm{(}\sqrt{\mathrm{1}\mathrm{+}\mathit{ϵ}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}{\mathit{ϵ}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\right]& & \end{array}$(B.2)where ${\mathit{\gamma }}_{\max }\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}{\mathit{\beta }}_{\mathrm{DGP}}}$ and $\begin{array}{ccc}\mathit{ϵ}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{8}\mathrm{(}{\mathit{r}}_{\mathrm{c}}{\mathit{H}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}}{\mathrm{9}{\mathit{\beta }}_{\mathrm{DGP}}^{\mathrm{2}}}\frac{{\mathrm{\Delta }}_{\mathrm{vir}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{c}\mathrm{\right)}}{{\mathit{x}}^{\mathrm{3}}{\mathit{g}}_{\mathrm{NFW}}\mathrm{\left(}\mathit{cx}\mathrm{\right)}}\mathrm{·}& & \end{array}$(B.3)

All Tables

All Figures

	Fig. 1 γ vs. R/R200 for NFW halos with masses ranging from 109 M⊙h-1 to 1014.5 M⊙h-1. The left panel shows the Symmetron solution using (assb,L,β) = (0.7, 1 Mpc h-1, 1), the central panel the Hu & Sawickif(R) solution with (| fR0 | ,n) = (10-6, 1), and the right panel the nDGP solution with rc = 1.
In the text

	Fig. 2 Semi-analytical solutions for γ for the Symmetron (left) and f(R) (right) using the same parameter values (and color coding) as in Fig. 1.
In the text

Fig. 3 vs halo mass M. The left panel shows the Symmetron curves with (assb,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 | fR0 | = (−4, −5, −6, −7) (in red, blue, green, and purple, respectively) and the right panel the nDGP solution with rc = (0.5, 1 ,3) (from top to bottom).

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.
In the text

	Fig. 5 versus M/μ200 for the Hu-Sawicki 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 fR0. μ200 for this model is given by Eq. (38).
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.
In the text

	Fig. 7 Semi-analytical predictions for versus M/μ200 for the Hu-Sawicki 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.
In the text

Fig. 8 γ(M) for NFW halos with masses ranging from 109 M⊙h-1 to 1015 M⊙h-1 using Eq. (52) with $\mathit{G}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{\left(}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{\zeta z}{\mathrm{)}}^{\mathit{n}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{\zeta z}}$. 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.

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 semi-transparent areas). The red line shows the space spanned by the Hu-Sawicki f(R) model when varying | fR0 | 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, assb = 1/2. The orange line indicates the region for nDGP models with rc ∈ [0.5, 10] (left to right) with the associated semi-transparent area spanning wider than the plot range. See Sect. 4.1 for details.

In the text