© ESO, 2014

1. Introduction

Although there is overwhelming evidence that dark energy exists, its nature is still a mystery. One possible explanation is that dark energy is some type of still unknown field. The other is that general relativity is not fully correct and therefore requires some modification which can account for the phenomena attributed to dark energy. Pursuing the latter approach, often titled the modification of gravity, intense efforts have been made to establish new theories (see e.g. Clifton et al. 2012, for a theoretical overview). The main problems concerning the development of new gravitational theories are theoretical (e.g. the existence of ghosts), as well as observational, meaning the constraints from local experiments are very tight. One explanation for the possible variations that have not been discovered are so-called screening mechanisms, where the local constraints are fulfilled by hiding the additional force in dense environments such as the solar system.

In this work, two different families of scalar-tensor theories with screening are discussed: the symmetron model (Hinterbichler & Khoury 2010) and an f(R)-gravity model that uses chameleon screening (Hu & Sawicki 2007). The first model works through the dependency of the vacuum expectation value of the scalar field on the local matter density, whereas in the latter, the mass of the scalar field is linked to the matter density. Both mechanisms lead to the same outcome. In low-density regions, a fifth force is mediated, and in high-density regions this force is suppressed recovering general relativity. This behaviour requires astrophysical tests in order to examine the existence of such scalar fields.

Due to the non-linearity of the screening, N-body simulations are a preferable tool for observational predictions. For both of these models, this has been done before with the main interest in large-scale imprints, which consist mostly in deviating the matter power spectra and halo mass functions. See Brax et al. (2012), Davis et al. (2012), Winther et al. (2012), and Llinares & Mota (2013) for the symmetron case and Li et al. (2012), Brax et al. (2013), Li & Hu (2011), and Zhao et al. (2011) for the chameleon. Simulations are also able to predict individual halo properties, such as density or velocity dispersion profiles, which are well established observables. These observables can be analysed using gravitational lensing (Oguri et al. 2012; Johnston & Sheldon 2007; Okabe et al. 2013). Possible constraints for the f(R) model have been established by Lombriser et al. (2012) and Schmidt (2010).

Another possible test of gravity on scales below 1 Mpc is the use of gravitational redshift, since the wavelength shift of light is directly proportional to the depth of the potential well of the clusters of galaxies. This has been observed for the SDSS survey data by Wojtak et al. (2011) who compared the data points to an analytical prediction for general relativity based on the NFW profile (Navarro et al. 1995), an analytical TeVeS prediction¹, and a semi-analytical profile from Schmidt (2010). For the last, the NFW prediction has been boosted by the factor 4/3, which is the maximum enhancement of gravitational strength according for f(R) gravity.

The first detection of gravitational redshift in clusters was done by Wojtak et al. (2011) using 7800 clusters from the SDSS data set. A few corrections to their predictions were made later on (Zhao et al. 2013; Kaiser 2013) and pointed out that owing to the relative motion of the galaxies, an additional component exists because of time dilation. This term is called the transverse Doppler effect, and it has opposite sign as v_g. This and other effects, e.g., a changed redshift as a result of relativistic beaming, have been analysed by Kaiser (2013), all of them should be considered when doing a proper analysis of observational data. Because the scope of this work is, however, the deviation of v_g in screened modified gravity models with respect to ΛCDM, these additional effects have not been studied and have not been included in this work.

In this work, we used data coming from a set of eight N-body simulations that were run with the Isis code covering a variety of model parameters and presented in Llinares et al. (2013). The box size and resolution were chosen such that it is possible to analyse properties of haloes that correspond to groups and clusters of galaxies. We study the gravitational redshift profiles in the objects found in the simulations in the three gravitational models and show that such comparisons and predictions have to be made with extreme caution when dealing with screening models. During the analysis, special emphasis was put on determining the virialisation state of the haloes by taking the energy of the scalar field into account.

The outline of the paper is as follows. In Sect. 2 we review the concept of scalar-tensor theories with screening and with particular focus on the symmetron and the Hu-Sawicky f(R)-model. In Sect. 3 we describe the simulations and the methods employed for extracting information from the data sets. In particular, the halo selection and virialisation are discussed. In Sects. 4 and 5 we present and discuss the results. Finally, we conclude in Sect. 6.

2. Models

This section gives a general description of scalar-tensor theories with screenings mechanisms. Furthermore, the particular models that are the subjects of this paper are briefly reviewed. The implementation of the equations in the N-body code that was used to run the simulations is described in detail in Llinares et al. (2013).

2.1. Scalar-tensor theories with screening

The action of a scalar field φ in a scalar-tensor theory is given by $\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{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\nabla }}_{\mathit{\mu }}\mathit{\phi }{\mathrm{\nabla }}^{\mathit{\mu }}\mathit{\phi }\mathrm{-}\mathit{V}\mathrm{\left(}\mathit{\phi }\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)}& & \end{array}$(1)where the matter fields ψ⁽ⁱ⁾ couple to the Jordan frame metric given by $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{g}}_{\mathit{\mu \nu }}\end{array}\mathrm{\equiv }{\mathit{A}}^{\mathrm{2}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}{\mathit{g}}_{\mathit{\mu \nu }}\mathit{.}& & \end{array}$(2)This leads to an equation of motion for the scalar field $\begin{array}{ccc}\Lambda \mathit{\phi }\mathrm{=}{\mathit{V}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{-}{\mathit{A}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathit{T}& & \end{array}$(3)where T is the trace of the energy-momentum tensor. The right hand-side of this equation can be identified with an effective potential: $\begin{array}{ccc}{\mathit{V}}_{\mathrm{eff}}\mathrm{=}\mathit{V}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathit{A}\mathrm{\right(}\mathit{\phi }\mathrm{)}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{T}\mathit{.}& & \end{array}$(4)As in a matter-dominated background T ≈  − ρ_m, the value of the minimum depends on the local matter density. The screening nature of a model is now given if V_eff has a minimum close to zero for high-density backgrounds.

To describe the behaviour of the field in the environment of matter, we use the Newtonian gauge metric given by $\begin{array}{ccc}\hspace{0.17em}\mathrm{d}{\mathit{s}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathrm{\Phi }\mathrm{)}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\mathrm{\Phi }\mathrm{)}\mathrm{(}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{2}}\mathrm{+}\hspace{0.17em}\mathrm{d}{\mathit{y}}^{\mathrm{2}}\mathrm{+}\hspace{0.17em}\mathrm{d}{\mathit{z}}^{\mathrm{2}}\mathrm{)}\mathit{.}& & \end{array}$(5)The geodesic equation in this particular coordinates becomes $\begin{array}{ccc}\begin{array}{c}{¨}\\ {x}\end{array}\mathrm{+}\mathrm{2}\mathit{H}{ẋ}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{\nabla }\left(\mathrm{\Phi }\mathrm{+}\log \mathit{A}\right)\mathrm{=}\mathrm{0}& & \end{array}$(6)where the quasi-static limit has been applied (i.e. the time derivatives of the field were assumed to be much smaller than its spatial variation). The last term of the equation can be interpreted as an effective force potential, and the arising fifth force is given by $\begin{array}{ccc}{\begin{array}{c}{¨}\\ {x}\end{array}}_{\mathrm{Fifth}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{\nabla }\log \mathit{A}\mathit{.}& & \end{array}$(7)

2.2. Symmetron

In the symmetron model (Hinterbichler & Khoury 2010; Hinterbichler et al. 2011; Olive & Pospelov 2008), the expectation value of the scalar field depends on the local matter density. In high-density regions it is zero, but when the matter density is lower than a given threshold, the symmetry of the potential is broken, and the potential acquires a minimum away from zero, leading to an additional force in these regions.

The requirement for the symmetron is that both A(φ) and V(φ) are symmetric under φ →  − φ, and therefore a commonly considered coupling is $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{{\mathit{\phi }}^{\mathrm{2}}}{\mathrm{2}{\mathit{M}}^{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\left(\frac{{\mathit{\phi }}^{\mathrm{4}}}{{\mathit{M}}^{\mathrm{4}}}\right)& & \end{array}$(8)where M is a mass scale. The simplest potential can be stated as $\begin{array}{ccc}\mathit{V}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\mu }}^{\mathrm{2}}{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}\mathit{\lambda }{\mathit{\phi }}^{\mathrm{4}}\mathit{,}& & \end{array}$(9)where μ is a mass scale and λ a dimensionless parameter. This definition leaves the effective potential (4) for non-relativistic matter as $\begin{array}{ccc}{\mathit{V}}_{\mathrm{eff}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathit{\rho }}{{\mathit{M}}^{\mathrm{2}}}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\right){\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}{\mathit{\phi }}^{\mathrm{4}}\mathit{,}& & \end{array}$(10)which has a minimum at φ = 0 for densities ρ > ρ_ssb ≡ μ²M² and two minima if ρ < ρ_ssb. Particularly, the two minima are ${\mathit{\phi }}_{\mathrm{0}}\mathrm{=}\mathrm{\pm }\mathit{\mu }\mathit{/}\sqrt{\mathit{\lambda }}$ in vacuum.

The geodesic equation for this model turns into $\begin{array}{ccc}\begin{array}{c}{¨}\\ {x}\end{array}\mathrm{+}\mathrm{2}\mathit{H}{ẋ}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{\nabla \Phi }\mathrm{+}\frac{\mathrm{1}}{\mathrm{(}\mathit{Ma}{\mathrm{)}}^{\mathrm{2}}}\mathit{\phi }\mathrm{\nabla }\mathit{\phi }\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(11)hence, ${\begin{array}{c}{}_{{¨}}\\ {x}\end{array}}_{\mathrm{Fifth}}\mathrm{\propto }\mathit{\phi }\mathrm{\nabla }\mathit{\phi }$. Consequently, in a high-density region, the fifth force is screened as required.

Similar to Winther et al. (2012), instead of the original parameters (μ,M,λ), more physical parameters were introduced that are all linked to the properties of the scalar field in vacuum. Firstly, the range of the field in vacuum $\begin{array}{ccc}\mathit{L}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}}\mathit{\mu }}\mathrm{;}& & \end{array}$(12)secondly, the expansion factor for which the symmetry is broken in the background level $\begin{array}{ccc}{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{m}\mathrm{0}}{\mathit{\rho }}_{\mathrm{c}\mathrm{0}}}{{\mathit{\mu }}^{\mathrm{2}}{\mathit{M}}^{\mathrm{2}}}\mathrm{;}& & \end{array}$(13)and finally a dimensionless coupling constant $\begin{array}{ccc}\mathit{\beta }\mathrm{=}\frac{\mathit{\mu }{\mathit{M}}_{\mathit{Pl}}}{{\mathit{M}}^{\mathrm{2}}\sqrt{\mathit{\lambda }}}\mathrm{·}& & \end{array}$(14)In addition, the scalar field itself is normalised with its vacuum expectation value: $\begin{array}{ccc}\mathit{\chi }\mathrm{\equiv }\mathit{\phi }\mathit{/}{\mathit{\phi }}_{\mathrm{0}}\mathrm{=}\frac{\mathit{\phi }{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}}{\mathrm{6}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{M}}_{\mathit{Pl}}{\mathit{L}}^{\mathrm{2}}\mathit{\beta }{\mathrm{\Omega }}_{\mathrm{m}\mathrm{0}}}\mathrm{·}& & \end{array}$(15)With these redefinitions, the geodesic equation can be rewritten to $\begin{array}{ccc}\begin{array}{c}{¨}\\ {x}\end{array}\mathrm{+}\mathrm{2}\mathit{H}{ẋ}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{\nabla \Phi }\mathrm{+}\frac{\mathrm{6}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{L}}^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}{\mathit{a}}_{\mathrm{ssb}}^{\mathrm{3}}}\mathit{\chi }\mathrm{\nabla }\mathit{\chi }\mathrm{=}\mathrm{0}\mathrm{·}& & \end{array}$(16)The equation of motion of the scalar field is in this case $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\chi }\mathrm{=}\frac{{\mathit{a}}^{\mathrm{2}}}{\mathrm{2}{\mathit{L}}^{\mathrm{2}}}\left[{\left(\frac{{\mathit{a}}_{\mathrm{ssb}}}{\mathit{a}}\right)}^{\mathrm{3}}\mathit{\eta \chi }\mathrm{-}\mathit{\chi }\mathrm{+}{\mathit{\chi }}^{\mathrm{3}}\right]& & \end{array}$(17)where η ≡ ρ_ma³/(Ω_m0ρ_c0) is the matter density in terms of the background density.

2.3. Hu-Sawicky f(R)-model

In the f(R) model by Hu & Sawicki (2007), the Jordan frame action is given by $\begin{array}{ccc}\mathit{S}& & \mathrm{=}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{4}}\mathit{x}\{\sqrt{\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{g}\end{array}}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{+}\mathit{f}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}{\mathrm{)}}^{\mathrm{)}}\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}$(18)where $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathit{m}}^{\mathrm{2}}\frac{{\mathit{c}}_{\mathrm{1}}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{/}{\mathit{m}}^{\mathrm{2}}{\mathrm{)}}^{\mathit{n}}}{\mathrm{1}\mathrm{+}{\mathit{c}}_{\mathrm{2}}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{/}{\mathit{m}}^{\mathrm{2}}{\mathrm{)}}^{\mathit{n}}}\mathrm{·}& & \end{array}$(19)Here, c₁, c₂, m, and n are positive constants and m is specified as ${\mathit{m}}^{\mathrm{2}}\mathrm{=}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{m}\mathrm{0}}$. Demanding ΛCDM background evolution (c₁/c₂m² = 2Λ) and taking the high curvature limit, which is an expansion around ${\mathit{m}}^{\mathrm{2}}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{\to }\mathrm{0}$, and introducing the current background curvature f_R0 ≡    df/$\hspace{0.17em}\mathrm{d}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}{|}_{\mathit{a}\mathrm{=}{\mathit{a}}_{\mathrm{0}}}$, Eq. (18) can be written as $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{16}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{-}\frac{{\mathit{f}}_{\mathit{R}\mathrm{0}}}{\mathit{n}}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}_{\mathrm{0}}^{\mathit{n}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\end{array}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{\mathit{n}}\end{array}}\mathrm{·}& & \end{array}$(20)The f(R) models can be treated as scalar-tensor theories, via a conformal transformation: $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{g}}_{\mathit{\mu \nu }}\end{array}\mathrm{=}{\mathrm{e}}^{\mathrm{2}\mathit{\beta \phi }\mathit{/}{\mathit{M}}_{\mathit{Pl}}}{\mathit{g}}_{\mathit{\mu \nu }} \mathrm{with}\mathit{\beta }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{6}}}\mathrm{·}& & \end{array}$(21)Under this transformation to the Einstein frame metric g_μν, a scalar field appears with $\begin{array}{ccc}{\mathit{f}}_{\mathit{R}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\beta \phi }\mathit{/}{\mathit{M}}_{\mathit{Pl}}}\mathrm{-}\mathrm{1}\mathrm{\approx }\mathrm{-}\frac{\mathrm{2}\mathit{\beta \phi }}{{\mathit{M}}_{\mathit{Pl}}}\mathit{,}& & \end{array}$(22)and its potential is given by $\begin{array}{ccc}\mathit{V}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{M}}_{\mathit{Pl}}^{\mathrm{2}}\frac{{\mathit{f}}_{\mathit{R}}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{-}\mathit{f}}{\mathrm{2}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{f}}_{\mathit{R}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(23)The effective linearised potential reads as $\begin{array}{ccc}{\mathit{V}}_{\mathrm{eff}}& \mathrm{=}\mathit{V}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{+}\frac{\mathit{\beta \phi }}{{\mathit{M}}_{\mathit{Pl}}}\mathit{\rho }\mathrm{·}& \end{array}$(24)By inserting the above definitions in Eq. (6), the geodesic equation for this model can be written as $\begin{array}{ccc}\begin{array}{c}{¨}\\ {x}\end{array}\mathrm{+}\mathrm{2}\mathit{H}{ẋ}\mathrm{+}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{\nabla }\left(\mathrm{\Phi }\mathrm{-}\frac{{\mathit{f}}_{\mathit{R}}}{\mathrm{2}}\right)\mathrm{=}\mathrm{0.}& & \end{array}$(25)Hence, the additional force in this case is ${\begin{array}{c}{}_{{¨}}\\ {x}\end{array}}_{\mathrm{Fifth}}\mathrm{\propto }\mathrm{\nabla }{\mathit{f}}_{\mathit{R}}\mathrm{\propto }\mathrm{\nabla }\mathit{\phi }$. Since the minimum of V_eff approaches zero with a growing ρ, the fifth force is screened in high-density regions as seen before.

The parameters of the Hu-Sawicky model discussed before, namely n and f_R0, were also used in the simulation. Basically, the effective potential and the conformal transformation described can be used directly to solve the modified geodesics and the scalar field evolution numerically as done for the symmetron. Sin the scalar field has a fixed sign in this case, it is customary to stabilise the numerical solution by forcing its sign to be unique (Oyaizu 2008). This can be made by including a new change of variables of the following form: $\begin{array}{ccc}{\mathrm{e}}^{\mathit{u}}\mathrm{\equiv }\mathrm{-}{\mathit{f}}_{\mathit{R}}{\mathit{a}}^{\mathrm{2}}\mathit{.}& & \end{array}$(26)This leads to the equation of motion for the f_R-field $\begin{array}{ccc}\mathrm{\nabla }\mathrm{·}\left(\exp \mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{\nabla }\mathit{u}\right)& \mathrm{=}& {\mathrm{\Omega }}_{\mathrm{m}}\mathit{a}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathrm{-}\mathrm{1}\mathrm{)}\\ & & \mathrm{-}{\mathrm{\Omega }}_{\mathrm{m}}\mathit{a}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}\mathrm{+}\mathrm{4}\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_{\mathrm{m}}}\right)\mathrm{\left(}\mathrm{\right|}{\mathit{f}}_{\mathit{R}\mathrm{0}}\mathrm{|}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{u}}{\mathit{n}\mathrm{+}\mathrm{1}}}\\ & & \mathrm{+}{\mathrm{\Omega }}_{\mathrm{m}}\mathit{a}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{a}}^{\mathrm{3}}\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_{\mathrm{m}}}\right)\mathrm{·}\end{array}$(27)

3. Simulations and analysis

3.1. The simulations

The simulations that we used for the analysis were run with the code Isis (Llinares et al. 2013), which is a modification of RAMSES (Teyssier 2002), which includes the development of the scalar fields. All simulation runs contained 512³ dark-matter-only particles with a mass of 9.26138 × 10⁹  M_⊙/h. The background cosmology is defined as (Ω_m0,Ω_Λ0,H₀) = (0.267,   0.733,   71.9). The simulation box has a side length at redshift zero of 256  h^-1 Mpc and the boundary conditions are periodic. The data sets are the snapshots taken at z = 0. All simulations were run with the same initial conditions, which were generated using Zeldovich approximation with standard gravity with the package Cosmics (Bertschinger 1995). In doing this, it was assumed that both extended models give fully screened fields before the initial redshift of the simulation. We refer the reader to Llinares et al. (2013) for further details on the simulations.

The parameters for the f(R) and the symmetron models are summarised in Table 1. In the case of the f(R)-gravity, the parameter n was fixed to one, while f_R0 took values from 10^-6, which resulted in hardly any deviation from ΛCDM, to 10^-4, which is on the border of violating cluster abundance constraints (Ferraro et al. 2011). For the symmetron model mainly the time of the symmetry breaking was varied while leaving the others constant except for one of the model that also has an increased coupling constant.

3.2. Halo selection

The halo identification was made using the halo finder Rockstar (Behroozi et al. 2013), a publicly available 6D FOF code. The boundary of the haloes was chosen to be R_200c (i.e. where the density falls below v = 200 times the critical density today). A proper definition of halo properties in screened modified gravity would require v to vary with the halo mass and the local environment, as spherical collapse analysis shows (Brax et al. 2010; Li & Efstathiou 2012). We found, however, that our results are stable against the change in v for the modified gravity models. For further analysis, we refer to this halo boundary as virial radius (i.e. R_v ≡ R_200c). Since we are dealing with gravitational effects, for which the total mass is the important quantity, we include the subhaloes of a given host as part of the main halo throughout the analysis. We defined the centre of the objects as the position of the particle that corresponds to the minimum gravitational potential. This choice is aimed at getting the difference in gravitational potential between the brightest cluster galaxy (BCG) and the rest of the cluster, thereby imitating an observational standpoint. This can, however, be a problem because the minimum of the gravitational potential is not necessarily the densest central region. To verify the findings, the results were also reproduced by taking the halo definition of the biggest subhalo from the Rockstar halo finder, and only minor deviations inside the virial radius were found.

The ΛCDM run contains around 73 000 dark matter haloes with more than 100 particles and circa 9200 with more than 1000 particles. The rest of the simulations have similar mass functions with small differences that come from the existence of the scalar field². An in-depth interpretation of the change in the halo mass functions in the symmetron and chameleon model was done by Brax et al. (2012, 2013) and is not the subject of this work.

3.3. Virialisation

Fig. 1 Distribution of the virialisation parameter βvir for haloes with mass M > 3 × 1013  M⊙   h-1. βvir using the potential energy obtained through direct summation WDS, which was defined in Eq. (37), is shown in orange with upward hatching, and βvir using W, defined in Eq. (33), is shown in blue with downward hatching. The vertical black lines visualise our criterion for virialisation, i.e. | βvir |  < 0.2.

To study halo properties, it is crucial not to mix haloes that are dynamically relaxed and those that have not reached such an equilibrium state yet because the two groups have a different distribution of halo properties (Shaw et al. 2006). Using the whole cluster sample and not separating these two groups may therefore lead to a skewed mean and higher noise. Especially when comparing the effect of different gravitational theories below Mpc scale – as done in this work – combining the two groups may lead to a false interpretation because it is not clear if the change happens on an actual halo property or if only the number of unrelaxed objects has changed.

By definition the energy portions of a virialised object fulfil the virial theorem: $\begin{array}{ccc}\mathrm{2}\mathit{T}\mathrm{+}\mathit{W}\mathrm{-}{\mathit{E}}_{\mathrm{S}}\mathrm{=}\mathrm{0}& & \end{array}$(28)with the potential energy W, the kinetic energy T, and the surface pressure term E_s. For a perfect fluid with density ρ and pressure p, these quantities are given by $\begin{array}{ccc}\mathit{W}\mathrm{=}{\mathrm{\int }}_{\mathit{V}}\mathrm{d}{}^{\mathrm{3}}{x}\hspace{0.17em}\mathit{\rho }{x}\mathrm{·}\begin{array}{c}{¨}\\ {x}\end{array}\mathit{,}& & \end{array}$(29)$\begin{array}{ccc}{\mathit{E}}_{\mathrm{S}}\mathrm{=}\underset{}{{\mathrm{\int }}^{\mathrm{S}}}\hspace{0.17em}\mathrm{d}{S}\mathrm{·}\mathrm{\left(}{x}\mathit{p}\mathrm{\right)}& & \end{array}$(30)and $\begin{array}{ccc}\mathit{T}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\underset{}{{\mathrm{\int }}^{\mathit{V}}}\mathrm{d}{}^{\mathrm{3}}{x}\hspace{0.17em}\mathit{p}\mathrm{·}& & \end{array}$(31)Since these relations are derived from the collisionless Boltzmann equation (e.g. Chandrasekhar 1981), they are true in any modified gravity model, hence also in the models we consider. Any modification in gravitational force enters through the total acceleration $\begin{array}{c}{{¨}}_{}\\ {x}\end{array}$, which is a sum of the Newtonian part and an acceleration induced by a fifth force.

The previous definitions must be discretised if one wants to apply them to N-body simulations. The kinetic energy term then yields $\begin{array}{ccc}\mathit{T}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}}{\mathit{m}}_{\mathit{i}}{{v}}_{\mathit{i}}^{\mathrm{2}}& & \end{array}$(32)where v_i is the relative velocity of particle i with respect to the halo velocity.

The potential energy was obtained using the acceleration of the particles: $\begin{array}{ccc}\mathit{W}\mathrm{=}\sum _{\mathit{i}}\frac{{\mathit{m}}_{\mathit{i}}}{\mathrm{2}}\mathrm{(}{{x}}_{\mathit{i}}\mathrm{-}{{x}}_{\mathrm{\mathscr{H}}}\mathrm{)}\mathrm{·}{\begin{array}{c}{¨}\\ {x}\end{array}}_{\mathit{i}}\mathit{,}& & \end{array}$(33)where the centre of the halo x_ℋ was defined as the minimum of the gravitational potential. Since the acceleration was calculated using the Newtonian and the fifth force, the definition ensures that the energy of the scalar field is not neglected. Both forces were obtained directly through the N-body code via an inverse cloud-in-cell (inverse CIC) smoothing, in the same way as done while running the simulations.

The E_S term was calculated following Shaw et al. (2006), who use an approximation that is motivated by the ideal gas law and reads as $\begin{array}{ccc}{\mathit{E}}_{\mathrm{s}}\mathrm{\approx }\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{90}}^{\mathrm{3}}{\mathit{p}}_{\mathrm{S}}\mathit{,}& & \end{array}$(34)where $\begin{array}{ccc}{\mathit{p}}_{\mathrm{S}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{{\sum }_{\mathit{i}}{\mathit{m}}_{\mathit{i}}{{v}}_{\mathit{i}}^{\mathrm{2}}}{\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi }\mathrm{(}{\mathit{R}}_{\mathrm{100}}^{\mathrm{3}}\mathrm{-}{\mathit{R}}_{\mathrm{80}}^{\mathrm{3}}\mathrm{)}}\mathrm{·}& & \end{array}$(35)Here, the sum in the nominator was carried out using the outermost 20% of the particles within R_vir, and the notation R_x is used to describe the x-percent quantile of this particle distribution.

Measuring the level of relaxation can be done using the virialisation parameter, which is defined as $\begin{array}{ccc}{\mathit{\beta }}_{\mathrm{vir}}\mathrm{\equiv }\frac{\mathrm{2}\mathit{T}\mathrm{-}{\mathit{E}}_{\mathrm{s}}}{\mathit{W}}\mathrm{+}\mathrm{1.}& & \end{array}$(36)We then defined a halo as being sufficiently relaxed if | β_vir |  < 0.2. This removes around half of all the haloes from each data sample. Mainly smaller haloes are removed. For the ΛCDM data, the three halo mass bins that will be adopted later on (i.e. log   (M   h/M_⊙) ∈ ((13,   13.5),  (13.5,   14),  (14,   14.5))), loose (19,  13,  18) percent of the haloes, which leaves (5060,  1544,  314) haloes. The three haloes in the ΛCDM data set with masses higher than 10^14.5  M_⊙   h^-1 are also unvirialised according to our definition. These rather large quantities of removed haloes are mainly due to our previously introduced halo definition after which subhaloes are merged into the host halo, so that the full phase-space analysis of Rockstar is partially lost.

Figure 1 shows the distribution of the virialisation parameter for haloes with masses M > 3 × 10¹³  M_⊙   h^-1. To illustrate that the energy of the scalar field must not be neglected, β_vir has been calculated not only with a potential energy as stated in Eq. (33) but also with a potential energy obtained through direct summation using only standard gravity: $\begin{array}{ccc}{\mathit{W}}_{\mathrm{DS}}\mathrm{=}\mathrm{-}\mathit{G}\sum _{\mathit{ij}}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{\mathrm{|}{{x}}_{{i}}\mathrm{-}{{x}}_{{j}}\mathrm{|}}\mathrm{·}& & \end{array}$(37)For ΛCDM, the fofr6 or the symm_A model, the deviation is not significant, but in the more extreme models (fofr4, symm_D) the direct summation method is clearly not sufficient. For them, the higher kinetic energy due to additional acceleration is not compensated through an equally increased potential energy and the β_vir-distribution is shifted to lower values.

3.4. Measuring the gravitational redshift from the simulations

In order to obtain the gravitational redshift, first the gravitational potential at each of the particles’ positions is needed. To achieve this, the N-body code was modified to interpolate several properties into the particles’ positions and output them. Amongst these properties were, for example, the matter density and the gravitational potential.

The measure of gravitational redshift was then defined as $\begin{array}{ccc}{\mathit{v}}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{\Delta \Phi }}{\mathit{c}}\mathrm{=}\frac{{\mathrm{\Phi }}_{\mathrm{c}}\mathrm{-}\mathrm{\Phi }}{\mathit{c}}\mathrm{·}& & \end{array}$(38)Here, Φ_c is the minimum gravitational potential. The particle that corresponds to this value was used (as in previous section) to define the centre of the halo. Therefore, v_g is always negative and can be interpreted as the gravitational blueshift of light seen by an observer located at the centre of the halo.

Instead of defining a single point as centre of the cluster and calculating the gravitational redshift with respect to that point, it is also possible to model the central galaxy as being spread out and take an average over gravitational potentials of most central particles as the reference point. This procedure is used by Kim & Croft (2004) and leads to a flattening of the profile. However, since they found out the effect is rather weak for dark-matter-only simulations, and our results are consistent with the NFW predictions we do not adopt the alternative procedure.

Our approach yields values for v_g(R_vir) of approximately − 7 km s^-1 and − 1.5 km s^-1 in the mass bins 10¹³ − 10^13.5  M_⊙/h and 10¹⁴ − 10^14.5  M_⊙/h, respectively. This is consistent with the findings of Wojtak et al. (2011) and Kim & Croft (2004), keeping in mind their altered halo mass bins.

Fig. 2 Additional acceleration $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Fifth}}$ due to the presence of a scalar field in the f(R) (top) and symmetron (bottom) models as a function of radius for three different mass bins. Also, the pure Newtonian gravitational force $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Newton}}$ is included from the ΛCDM data set (black, unfilled circles) as a comparison. The vertical line approximately corresponds to two grid cells in the finest refinement level. See text for details.

4. Results

The aim of the paper is to study possible differences in the gravitational redshift signal owing to the presence of a scalar field. However, since null geodesics are invariant under conformal transformations and, thus, identical in the Einstein and Jordan frame, the photons are not affected by the presence of scalar fields. Consequently, the gravitational redshift is not influenced directly by the scalar field, but maps out the change in matter clustering due to the presence of the fifth force.

To understand the change in the gravitational redshift profiles, we also analysed the magnitude of the fifth force. For each of the quantities under study, we obtained an estimation of the error by dividing the simulation box into eight sub-boxes and using the variance of the relative deviation calculated in each subset to fix the error of the mean.

4.1. Distribution of fifth force in the dark matter haloes

Figure 2 shows the absolute value of the acceleration $\mathrm{\left|}{\begin{array}{c}{}_{{¨}}\\ {x}\end{array}}_{\mathrm{Fifth}}\mathrm{\right|}$ defined by Eq. (7) as a function of radius at redshift z = 0. The forces and scalar field that are needed to evaluate $\mathrm{\left|}{\begin{array}{c}{}_{{¨}}\\ {x}\end{array}}_{\mathrm{Fifth}}\mathrm{\right|}$ were calculated in the same way as we did when calculating the virialisation state (i.e. using the same smoothing kernel that the N-body code Isis used during the simulations to interpolate quantities from the grid into the particle’s position). The adaptive mesh refinement of the code was used during these calculations, and so we reach the same resolution as during the simulations. The upper and bottom rows correspond to f(R) and symmetron simulations, respectively. We show results for three different mass bins for each model. The vertical line to the left of the panels corresponds to the resolution limit that was estimated as twice the size of the grid of the deepest refinement level normalised with the mean virial radius of each mass bin.

All the curves show a characteristic maximum, whose position is highly dependent on the model parameters and the mass of the haloes. In the f(R) case, the maximum of the fifth force profile moves towards larger radii when increasing mass and decreasing the only free parameter | f_R0 |. This is a direct consequence of the screening that is activated in a larger part of the haloes when decreasing | f_R0 |. For the symmetron model, an increase of β or z_SSB leads to a stronger fifth force. A greater β value increases the $\mathrm{\left|}\begin{array}{c}{}_{{¨}}\\ {{x}}_{\mathrm{Fifth}}\end{array}\mathrm{\right|}$ values by a constant factor (symm_A versus symm_C) while altering z_ssb changes the shape of the fifth force profile in general.

To better understand the mass dependence of the distributions, we show in Fig. 3 the absolute value of the acceleration $\mathrm{\left|}{\begin{array}{c}{}_{{¨}}\\ {x}\end{array}}_{\mathrm{Fifth}}\mathrm{\right|}$ as a function of halo mass as measured at two different radii. The information presented in the plots was obtained from two spherical shells of radius 0.1  R_v and R_v and a thickness of 0.02 R_v. The upper panels of the figure show the results from the f(R) simulations, while the bottom panels correspond to the symmetron data. From this figure, it is clear that different sized haloes are variably affected by the fifth force. The model parameters lead to a characteristic halo mass range beyond which the fifth force is screened. Below this range, the haloes are not enough dense for the fifth force to be activated. In the case of the f(R) model, we find that the low mass end of the force distribution is completely insensitive to changes in the only free parameter f_R0. On the contrary, we find a strong dependence in the high mass end of the distribution: the larger f_R0, the higher the masses that are screened and thus affected by the fifth force. This produces a shift in the maximum of the distributions towards higher masses when f_R0 is increased.

Fig. 3 Additional acceleration $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Fifth}}$ due to the presence of a scalar field in the f(R) (top) and symmetron (bottom) models as a function of halo mass for two different radial bins. As in Fig. 2, the black unfilled circles represent $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Newton}}$ as a comparison.

Fig. 4 Relative deviation in gravitational redshift between the scalar-tensor models and ΛCDM as a function of radius. The upper panels correspond to the f(R) model and the lower to the symmetron. The vertical line approximately corresponds to two grid cells in the finest refinement level.

The dependence of the symmetron fifth force on the redshift of symmetry breaking z_SSB is similar to the dependence of this quantity on f_R0 in the Hu-Sawicky model: the higher the symmetry breaking, the higher the masses that are unscreened and affected by the fifth force. Since the low mass end of the distribution is also rather insensitive to changes in z_SSB, there is a displacement of the maximum of the distributions towards high masses when increasing z_SSB. By comparing simulations symm_A and symm_C, we can test the dependence of the forces when changing the coupling constant β. We find a rather insensitive behaviour. The only changes are in the normalisation and are given by the dependence with β in Eq. (16).

Fig. 5 Relative deviation in gravitational redshift between the scalar-tensor models and ΛCDM as a function of halo mass for two different radial bins. The upper panels correspond to the f(R) model and the lower to the symmetron.

4.2. Gravitational redshift

The aim of the paper has been to study gravitational redshift. As discussed in the introduction, in the family of alternative theories that we are treating in this paper, the scalar field does not affect the energy of photons, and thus, any difference in the gravitational redshift prediction will come from differences in the matter distribution. Naively, one would expect that the presence of a fifth force will increase the clustering and, thus, produce deeper potential wells, with the consequence of an increase in the gravitational redshift at all masses. Our simulations show that the situation is a bit more complex.

Figure 4 shows the relative deviation in the gravitational redshift v_g with respect to ΛCDM as a function of radius. The upper and lower panels show the f(R) and symmetron results, respectively. As previously in Fig. 3, the vertical black line denotes and estimation of the resolution limit. Since the fifth force causes additional clustering, which affects the gravitational redshift profile of the clusters, we expect to find similar dependencies on the model parameters to those in the previous section. This clearly applies to the f(R) results where we find, as expected, that the smaller | f_R0 |, the smaller the deviation from ΛCDM. Especially in the | f_R0 |  = 10^-6 case, the result is basically indistinguishable from ΛCDM. On the other hand, for the more extreme cases of | f_R0 |  = 10^-4 and | f_R0 |  = 10^-5, the relative deviation is much greater, in particular for small radii. In the case of the symmetron model, the imprint in the variation of the gravitational redshift profiles is not as strong as in the f(R) results, which is a consequence of the fact that the symmetron is a more effective screening mechanism. The maximum deviation in this case is around 15% for the symm_D model for haloes with masses between 10¹⁴ and 10^14.5  M_⊙   h^-1. As expected, this deviation decreases for later symmetry breaking times. In particular, the symm_A model with z_ssb = 1 is basically indistinguishable from ΛCDM. For most of the parameter sets and halo mass bins analysed, the results show a stronger gravitational redshift profile in the modified gravity models, following the expected behaviour. However, the models symm_B, symm_D and fofr5 show a signal that contradicts the naive expectations (i.e. show a negative correction towards a less important gravitational redshift). Before discussing the reason for this happening, we briefly discuss the dependence of the gravitational redshift on mass.

The three halo mass bins displayed in the panels of Fig. 4 already suggest that the amount of deviation in gravitational redshift depends on the cluster mass. For instance, the fofr5 model shows the largest enhancement in gravitational redshift in the low mass bin, while the fofr4 becomes dominant for the highest mass haloes. This mass dependency of the deviation in gravitational redshift is shown more clearly in Fig. 5, where we present the relative deviation with respect to Einstein gravity as a function of mass and for two different radii. Here, we find mass ranges for which the deviation is greatest. As previously noticed for the characteristic mass ranges noticed in Fig. 3, the position of the peak depends on the model parameters. For example, in the fofr6 data the maximum deviation is for haloes with masses around 2 × 10¹²  M_⊙/h, and for fofr5 it is around 3 × 10¹³M_⊙/h. The fofr5 model does not present a peak, but a continuous growth towards the highest mass haloes. Larger box sizes are needed to map the larger clusters and to confirm that the signal is also the same for this model. Similarly, the symmetron curves also possess maxima whose positions correspond to the ones of the peaks in Fig. 3.

5. Discussion

It is well known that scalar tensor models increase the clustering rate (see for instance results from N-body simulations in Brax et al. 2013, 2012; Li et al. 2012). In particular, we find:

• f(R) results: the lower the value of | f_R0 |, the smaller are the affected haloes, and the fifth force gets reduced. Since f_R0 controls directly, the deviation from ΛCDM in the action (see Eq. (20)) and the range of the field in vacuum is λ_φ   =   1/$\sqrt{{\mathit{V}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}}\mathrm{\propto }\sqrt{\mathrm{\left|}{\mathit{f}}_{\mathit{R}\mathrm{0}}\mathrm{\right|}}$, which is the expected result.

• symmetron results: as expected from Eq. (16) a lower value of a_ssb or a higher value of β results in a stronger fifth force and a generally greater deviation from the ΛCDM cluster profiles. In this case the range of the field in vacuum L was not altered, and consequently, the size of the affected haloes did not change dramatically. Still, we find lowering a_ssb leads to bigger haloes being affected by the fifth force, since the symmetron force has more time to influence the matter clustering.

The most common expectation about the change in the gravitational redshift through this additional clustering is to observe an enhanced signal. This idea is motivated by the fact that additional clustering should lead to deeper potential wells, hence to a stronger gravitational redshift signal. Mostly, this is the major imprint observed in the f(R) models studied, and the same behaviour was expected in the symmetron models. However, depending on the halo mass and model parameters, the matter density can be enhanced in the outskirts of the halo. Keeping in mind the halo mass is effectively fixed in one mass bin, this result is not surprising as missing mass in one region has to be compensated by additional mass in another region. Such an altered halo density profile means the potential well becomes shallower in the central regions and a deeper in the outskirts. Consequently, the resulting deviation in gravitational redshift is negative. Still, this connection between clustering in the central regions and stronger gravitational redshift or clustering in the outer regions and a negative deviation does not explain the opposite imprint obtained for various parameter sets. Therefore, the question we tried to answer is what mechanism controls whether the clustering happens in the central region or in the outskirts of a dark matter halo?

To answer this question, the position of the additional clustering is analysed. This position is related to the radius of the maximal fifth force R_max (see Fig. 2) the following way. For R ≳ R_max, where the fifth force decreases because a lower matter density leads to a less affected φ-value, the additional mass-flux towards the centre is greater than the supplementary inflow of mass towards R_max. For R ≲ R_max, where the fifth force decreases due to a flattening of the halo density profile or because of screening, the opposite is the case. This means that the additional clustering happens around R_max. The extent of the additional clustering depends on the overall strength of the fifth force and the steepness of the slopes around R_max. In summary, an R_max close to the resolution limit leads to a deepening of the potential well and therefore to a positive deviation in gravitational redshift, whereas a greater value of R_max indicates an additional clustering in the outskirts of the halo, hence a negative deviation with respect to ΛCDM. Another way of phrasing it is that each set of model parameters results in a “preferred density”, that is, a matter density for which the fifth force is maximal.

This reasoning can most easily be tested by observing the dependency of the fifth force and the deviation of the gravitational redshift for the symm_D model since R_max ≳ 0.1  R_v in this case. Here, clustering the centre is effectively weakened for clusters with mass ≳10¹³  M_⊙   h^-1. Consequently, gravitational redshift is decreased in these haloes. However, the explanation is supported by several additional features in the results. For instance, by decreasing z_ssb, the effect should still persist in a weaker form – which it does in the symm_B data set. A special position is taken by the fofr4 data: the increase in the number of massive clusters is immense, suggesting the same or even stronger merging of clusters. On the other hand, the force enhancement close to the centre has not yet reached its maximum, even for the largest haloes of the data sample (top left panel of Fig. 3) opposed to the symm_D case, where this maximum lies around 10¹³  M_⊙   h^-1 for R ≈ 0.1  R_v. This means that, in the fofr4 case, the “effective preferred density” is even higher than the one present in the centre of the biggest clusters. Therefore, the clustering is dominant in the central parts, leading to a higher gravitational redshift as seen in Fig. 4. It is more appropriate to compare the symm_D with the fofr5 results owing to their greater similarity in the fifth force profiles. The fifth acceleration in the central region is, however, more present for haloes in the mass range 10^13.5 − 10¹⁴  M_⊙   h^-1 in the fofr5 data. For the next bigger haloes, ${\begin{array}{c}{{¨}}_{}\\ {x}\end{array}}_{\mathrm{Fifth}}$ is negligible in the central regions in both models. This explains the resulting positive deviation in gravitational redshift for the first halo-mass bin mentioned and the negative deviation in the second (Fig. 4). Of course, the different other characteristics of the symmetron and the chameleon model (e.g. the time dependence of the fifth force) have to be taken into account if the data is to be analysed in more detail. The explanation provided can, after all, clarify the partially opposite results obtained.

6. Conclusions

We performed a set of eight high-resolution simulation runs with 512³ particles and a box length of 256 Mpc h^-1 each. The data was used to analyse the density and gravitational redshift profiles of virialised clusters. For this purpose, the virialisation parameter was first computed for each halo with the inclusion of the scalar field energy. Then, in order to calculate the gravitational redshift of all particles, the value of the gravitational potential at the particle’s position was subtracted from the halo’s minimal gravitational potential. Ultimately, the modified gravity results could be compared to the ΛCDM values. These steps were not only undertaken for the gravitational redshift and fifth force but also for the matter density, particle density, and velocity dispersion (not shown). These quantities allowed us to analyse and explain the gravitational redshift results. The overall finding from the numerical results is that in both analysed models, the deviation from ΛCDM can vary enormously depending on the choice of parameters and the analysed halo size.

The main goal of this work was to study gravitational redshift profiles of virialised haloes in both the symmetron model and chameleon f(R)-gravity using N-body simulations. To obtain the state of virialisation, several methods were analysed, and the most adequate one – using the total acceleration of the particles – was employed. This allowed us to calculate the virialisation parameter taking the energy of the scalar field into account. We note, however, that the qualitative results can be reproduced when including all haloes. This leads merely to an increase in noise, especially in the halo outskirts. In spite of this, we chose to only include the virialised haloes in our analysis to give a more conservative prediction.

We found the results to be highly dependent on the halo mass, which means the consideration of multiple mass ranges is crucial when analysing observational data. In particular, three possible regimes were identified:

• (i)

  No deviation from ΛCDM. If the halo is fullyscreened (either through self-screening or environmentalscreening) or is too small to affect the scalar field, no deviation inthe density and, consequently, in the gravitational redshiftprofiles was observed.

• (ii)

  Enhanced gravitational redshift. This is the expected result since the fifth force leads to additional clustering in the centre of the halo. Additional matter in the central region means a deeper potential well and, therefore, a positive deviation compared to ΛCDM.

• (iii)

  Weaker gravitational redshift. If the fifth force is screened in the central region of the halo but not in the outskirts, a negative deviation can appear. The additional force leads to a matter overdensity in the outer regions while – with respect to the halo mass – the inner regions are less dense than in the ΛCDM case.

This shows that simply assuming the change in the gravitational redshift is equal to the maximum possible change in the gravitational constant, as frequently done when predictions are made, is not sufficient. Instead, the prediction obtained by N-body simulations as in this work should be used.

Croft (2013) shows error predictions for the full SDSS/BOSS (Aihara et al. 2011), BigBOSS (Schlegel et al. 2009), and Euclid (Laureijs et al. 2011) data sets. Accordingly, BigBOSS and Euclid should be able to map out the amplitude of the v_g curve with 6.5% and 4% precision, respectively. These values were obtained using a different technique than our approach. Instead of taking the difference between the gravitational redshift of the central galaxy (BCG) and each particle in a cluster, the data was compared pairwise. As a result, they cannot be transferred directly. In addition, the binning of data in halo mass, which is necessary to detect signatures due to screened gravity, will naturally increase the error. Nevertheless, the planned new sky surveys will allow restriction of the parameter space of some modified gravity models.

Acknowledgments

M.B.G. would like to thank Michael Kopp for his useful comments on the manuscript. C.L. and D.F.M. thank the Research Council of Norway FRINAT grant 197251/V30. D.F.M. is also partially supported by project CERN/FP/123618/2011 and CERN/FP/123615/2011. The simulations were performed on the NOTUR Clusters HEXAGON and STALLO, the computing facilities at the University of Bergen and Tromsø in Norway.

References

All Tables

All Figures

Fig. 1 Distribution of the virialisation parameter βvir for haloes with mass M > 3 × 1013  M⊙   h-1. βvir using the potential energy obtained through direct summation WDS, which was defined in Eq. (37), is shown in orange with upward hatching, and βvir using W, defined in Eq. (33), is shown in blue with downward hatching. The vertical black lines visualise our criterion for virialisation, i.e. | βvir |  < 0.2.

In the text

Fig. 2 Additional acceleration $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Fifth}}$ due to the presence of a scalar field in the f(R) (top) and symmetron (bottom) models as a function of radius for three different mass bins. Also, the pure Newtonian gravitational force $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Newton}}$ is included from the ΛCDM data set (black, unfilled circles) as a comparison. The vertical line approximately corresponds to two grid cells in the finest refinement level. See text for details.

In the text

Fig. 3 Additional acceleration $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Fifth}}$ due to the presence of a scalar field in the f(R) (top) and symmetron (bottom) models as a function of halo mass for two different radial bins. As in Fig. 2, the black unfilled circles represent $\mathrm{|}\begin{array}{c}{}_{{¨}}\\ {x}\end{array}{\mathrm{|}}_{\mathrm{Newton}}$ as a comparison.

In the text

	Fig. 4 Relative deviation in gravitational redshift between the scalar-tensor models and ΛCDM as a function of radius. The upper panels correspond to the f(R) model and the lower to the symmetron. The vertical line approximately corresponds to two grid cells in the finest refinement level.
In the text

	Fig. 5 Relative deviation in gravitational redshift between the scalar-tensor models and ΛCDM as a function of halo mass for two different radial bins. The upper panels correspond to the f(R) model and the lower to the symmetron.
In the text