Issue 
A&A
Volume 562, February 2014



Article Number  A9  
Number of page(s)  10  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201322403  
Published online  30 January 2014 
Gravitational redshift profiles in the f(R) and symmetron models
Institute of Theoretical Astrophysics, University of Oslo, Postboks 1029, 0315 Oslo, Norway
email: maxbg@astro.uio.no
Received: 30 July 2013
Accepted: 15 November 2013
Aims. We investigate the gravitational redshift in clusters of galaxies in the symmetron and HuSawicky f(R) models. The characteristic feature of both models is the screening mechanism that hides the fifth force in dense environments recovering general relativity.
Methods. We use Nbody simulations that were run with the code Isis, which includes scalar fields, to analyse the deviation of observables in modified gravity models with respect to ΛCDM.
Results. We find that the presence of the screening makes the deviation highly dependent on the halo mass. For instance, the f(R) parameters  f_{R0}  = 10^{5}, n = 1 cause an enhancement of the gravitational signal by up to 50% for haloes with masses between 10^{13} M_{⊙} h^{1} and 10^{14} M_{⊙} h^{1}. The characteristic mass range where the fifth force is most active varies with the model parameters. The usual assumption is that the presence of a fifth force leads to a deeper potential well and thus a stronger gravitational redshift. However, we find that in cases in which only the central regions of the haloes are screened, there could also be a weaker gravitational redshift.
Key words: gravitation / dark energy / galaxies: clusters: general / largescale structure of Universe / galaxies: halos / methods: numerical
© 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 socalled 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 scalartensor 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 lowdensity regions, a fifth force is mediated, and in highdensity 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 nonlinearity of the screening, Nbody simulations are a preferable tool for observational predictions. For both of these models, this has been done before with the main interest in largescale 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^{1}, and a semianalytical 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 Nbody 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 scalartensor theories with screening and with particular focus on the symmetron and the HuSawicky 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 scalartensor 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 Nbody code that was used to run the simulations is described in detail in Llinares et al. (2013).
2.1. Scalartensor theories with screening
The action of a scalar field φ in a scalartensor theory is given by (1)where the matter fields ψ^{(i)} couple to the Jordan frame metric given by (2)This leads to an equation of motion for the scalar field (3)where T is the trace of the energymomentum tensor. The right handside of this equation can be identified with an effective potential: (4)As in a matterdominated 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 highdensity backgrounds.
To describe the behaviour of the field in the environment of matter, we use the Newtonian gauge metric given by (5)The geodesic equation in this particular coordinates becomes (6)where the quasistatic 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 (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 highdensity 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 (8)where M is a mass scale. The simplest potential can be stated as (9)where μ is a mass scale and λ a dimensionless parameter. This definition leaves the effective potential (4) for nonrelativistic matter as (10)which has a minimum at φ = 0 for densities ρ > ρ_{ssb} ≡ μ^{2}M^{2} and two minima if ρ < ρ_{ssb}. Particularly, the two minima are in vacuum.
The geodesic equation for this model turns into (11)hence, . Consequently, in a highdensity 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 (12)secondly, the expansion factor for which the symmetry is broken in the background level (13)and finally a dimensionless coupling constant (14)In addition, the scalar field itself is normalised with its vacuum expectation value: (15)With these redefinitions, the geodesic equation can be rewritten to (16)The equation of motion of the scalar field is in this case (17)where η ≡ ρ_{m}a^{3}/(Ω_{m0}ρ_{c0}) is the matter density in terms of the background density.
2.3. HuSawicky f(R)model
In the f(R) model by Hu & Sawicki (2007), the Jordan frame action is given by (18)where (19)Here, c_{1}, c_{2}, m, and n are positive constants and m is specified as . Demanding ΛCDM background evolution (c_{1}/c_{2}m^{2} = 2Λ) and taking the high curvature limit, which is an expansion around , and introducing the current background curvature f_{R0} ≡ df/, Eq. (18) can be written as (20)The f(R) models can be treated as scalartensor theories, via a conformal transformation: (21)Under this transformation to the Einstein frame metric g_{μν}, a scalar field appears with (22)and its potential is given by (23)The effective linearised potential reads as (24)By inserting the above definitions in Eq. (6), the geodesic equation for this model can be written as (25)Hence, the additional force in this case is . Since the minimum of V_{eff} approaches zero with a growing ρ, the fifth force is screened in highdensity regions as seen before.
The parameters of the HuSawicky 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: (26)This leads to the equation of motion for the f_{R}field (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^{3} darkmatteronly particles with a mass of 9.26138 × 10^{9} M_{⊙}/h. The background cosmology is defined as (Ω_{m0},Ω_{Λ0},H_{0}) = (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.
Model parameters of the different simulation runs.
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^{2}. An indepth 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 × 10^{13} M_{⊙} h^{1}. β_{vir} using the potential energy obtained through direct summation W_{DS}, 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: (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 (29)(30)and (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 , 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 Nbody simulations. The kinetic energy term then yields (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: (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 Nbody code via an inverse cloudincell (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 (34)where (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 xpercent quantile of this particle distribution.
Measuring the level of relaxation can be done using the virialisation parameter, which is defined as (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 phasespace analysis of Rockstar is partially lost.
Figure 1 shows the distribution of the virialisation parameter for haloes with masses M > 3 × 10^{13} 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: (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 Nbody 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 (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 darkmatteronly 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^{13} − 10^{13.5} M_{⊙}/h and 10^{14} − 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 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 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 subboxes 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 defined by Eq. (7) as a function of radius at redshift z = 0. The forces and scalar field that are needed to evaluate were calculated in the same way as we did when calculating the virialisation state (i.e. using the same smoothing kernel that the Nbody 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 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 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 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 as a comparison. 
Fig. 4
Relative deviation in gravitational redshift between the scalartensor 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 HuSawicky 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 scalartensor 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^{14} 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^{12} M_{⊙}/h, and for fofr5 it is around 3 × 10^{13}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 Nbody 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/, 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 massflux 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^{13} 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^{13} 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^{14} M_{⊙} h^{1} in the fofr5 data. For the next bigger haloes, is negligible in the central regions in both models. This explains the resulting positive deviation in gravitational redshift for the first halomass 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 highresolution simulation runs with 512^{3} 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 Nbody 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 selfscreening 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 Nbody 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.
The analytical TeVeS profile used was later shown to be based on inappropriate assumptions (Bekenstein & Sanders 2012).
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
 Aihara, H., Allen de Prieto, C., An, D., et al. 2011, ApJS, 193, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Behroozi, P. S., Wechsler, R. H., & Wu, H.Y. 2013, ApJ, 762, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Bekenstein, J. D., & Sanders, R. H. 2012, MNRAS, 421, L59 [NASA ADS] [CrossRef] [Google Scholar]
 Bertschinger, E. 1995, unpublished [arXiv:astroph/9506070] [Google Scholar]
 Brax, P., Rosenfeld, R., & Steer, D. 2010, J. Cosmol. Astropart. Phys. 2010, 033 [Google Scholar]
 Brax, P., Davis, A.C., Li, B., Winther, H. A., & Zhao, G.B. 2012, JCAP, 10, 002 [Google Scholar]
 Brax, P., Davis, A., Li, B., Winther, H., & Zhao, G. 2013, JCAP, 04, 029 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1981, Hydrodynamic and hydromagnetic stability (Dover Publications, Incorporated) [Google Scholar]
 Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Croft, R. A. C. 2013, MNRAS, 434, 3008 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, A., Li, B., Mota, D., & Winther, H. A. 2012, ApJ, 748, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Ferraro, S., Schmidt, F., & Hu, W. 2011, Phys. Rev. D, 83, 063503 [NASA ADS] [CrossRef] [Google Scholar]
 Hinterbichler, K., & Khoury, J. 2010, Phys. Rev. Lett., 104, 231301 [NASA ADS] [CrossRef] [Google Scholar]
 Hinterbichler, K., Khoury, J., Levy, A., & Matas, A. 2011, Phys. Rev. D, 84, 103521 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004 [CrossRef] [Google Scholar]
 Johnston, D., Sheldon, E., Wechsler, R., et al. 2007, unpublished [arXiv:0709.1159] [Google Scholar]
 Kaiser, N. 2013, MNRAS, 435, 1278 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, Y.R., & Croft, R. A. C. 2004, ApJ, 607, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011 [arXiv:1110.3193] [Google Scholar]
 Li, B., & Efstathiou, G. 2012, MNRAS, 421, 1431 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Hellwing, W. A., Koyama, K., et al. 2012, MNRAS, 428, 743 [Google Scholar]
 Li, Y., & Hu, W. 2011, Phys. Rev. D, 84, 084033 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., & Mota, D. F. 2013, Phys. Rev. Lett., 110, 161101 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., Mota, D. F., & Winther, H. A. 2013 [arXiv:1307.6748] [Google Scholar]
 Lombriser, L., Schmidt, F., Baldauf, T., et al. 2012, Phys. Rev. D, 85, 102001 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1995, MNRAS, 275, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Oguri, M., Bayliss, M. B., Dahle, H. k., et al. 2012, MNRAS, 420, 3213 [NASA ADS] [CrossRef] [Google Scholar]
 Okabe, N., Smith, G., & Umetsu, K. 2013, ApJ, 769, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Olive, K. A., & Pospelov, M. 2008, Phys. Rev. D, 77, 043524 [NASA ADS] [CrossRef] [Google Scholar]
 Oyaizu, H. 2008, Phys. Rev. D, 78, 123523 [NASA ADS] [CrossRef] [Google Scholar]
 Schlegel, D. J., Bebek, C., Heetderks, H., et al. 2009 [arXiv:0904.0468] [Google Scholar]
 Schmidt, F. 2010, Phys. Rev. D, 81, 103002 [NASA ADS] [CrossRef] [Google Scholar]
 Shaw, L. D., Weller, J., Ostriker, J. P., & Bode, P. 2006, ApJ, 646, 815 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 364, 337 [Google Scholar]
 Winther, H. A., Mota, D. F., & Li, B. 2012, ApJ, 756, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Wojtak, R., Hansen, S. H., & Hjorth, J. 2011, Nature, 477, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, G.B., Li, B., & Koyama, K. 2011, Phys. Rev. D, 83, 044007 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, H., Peacock, J., & Li, B. 2013, Phys. Rev. D, 88, 043013 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1
Distribution of the virialisation parameter β_{vir} for haloes with mass M > 3 × 10^{13} M_{⊙} h^{1}. β_{vir} using the potential energy obtained through direct summation W_{DS}, 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 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 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 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 as a comparison. 

In the text 
Fig. 4
Relative deviation in gravitational redshift between the scalartensor 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 scalartensor 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.