A general framework to test gravity using galaxy clusters III: Observable-mass scaling relations in $f(R)$ gravity

We test two methods, including one that is newly proposed in this work, for correcting for the effects of chameleon $f(R)$ gravity on the scaling relations between the galaxy cluster mass and four observable proxies. Using the first suite of cosmological simulations that simultaneously incorporate both full physics of galaxy formation and Hu-Sawicki $f(R)$ gravity, we find that these rescaling methods work with a very high accuracy for the gas temperature, the Compton $Y$-parameter of the Sunyaev-Zel'dovich (SZ) effect and the X-ray analogue of the $Y$-parameter. This allows the scaling relations in $f(R)$ gravity to be mapped to their $\Lambda$CDM counterparts to within a few percent. We confirm that a simple analytical tanh formula for the ratio between the dynamical and true masses of haloes in chameleon $f(R)$ gravity, proposed and calibrated using dark-matter-only simulations in a previous work, works equally well for haloes identified in simulations with two very different -- full-physics and non-radiative -- baryonic models. The mappings of scaling relations can be computed using this tanh formula, which depends on the halo mass, redshift and size of the background scalar field, also at a very good accuracy. The results of this work will form part of a general framework for unbiased and self-consistent tests of gravity using data from present and upcoming galaxy cluster surveys. The new results can be used for accurate determination of the cluster mass using SZ and X-ray observables, which is vital for tests that use probes including cluster number counts.


INTRODUCTION
Clusters of galaxies are the largest gravitationally-bound objects in the Universe and, as tracers of the highest peaks of the primordial density perturbations, they are frequently used to test cosmological models. For example, global properties, such as the abundance of clusters, are expected to be highly sensitive to the values of cosmological parameters and to the strength of gravity on large scales, and they can also be predicted accurately with numerical simulations. These simulations can be combined with state-of-the-art observational data to place tight constraints on a wide range of cosmological models, including modifications to the theory of General Relativity (GR). The wealth of quality data that will be made available from ongoing and upcoming galaxy cluster surveys (e.g., Lawrence et al. 2007;Levi et al. 2013;Laureĳs et al. 2011;Ivezic et al. 2008;Jansen et al. 2001;Weisskopf et al. 2000;Merloni et al. 2012;Ade et al. 2016;Hasselfield et al. 2013) has the potential to revolutionise the precision of these constraints.
Accurate measurement of the cluster mass is vital in tests that make use of probes such as cluster number counts and the gas mass fraction. However, direct estimation of the cluster mass is gener-★ E-mail: m.a.mitchell@durham.ac.uk ally an expensive task which typically requires long exposure times, making it challenging to repeat for very large samples of clusters. A common alternative technique is to infer the cluster mass by using its relation to more easily measured observables including the X-ray luminosity and the integrated Sunyaev-Zel'dovich (SZ) flux. These mass 'proxies' are intrinsically linked to the gravitational potential and, as a result, they typically have a power-law mapping with the mass. A lot of work has therefore been carried out to study and calibrate these cluster mass scaling relations using various approaches: simulations that employ a comprehensive galaxy formation model (full physics), including cooling, star formation and feedback, have been used (e.g., Nagai et al. 2007;Fabjan et al. 2011;Truong et al. 2018;Cui et al. 2018); internal calibration has been employed via a joint likelihood analysis (e.g., Mantz et al. 2010Mantz et al. , 2014; selfcalibration can be achieved using additional observables, such as the clustering of clusters (e.g., Schuecker et al. 2003;Majumdar & Mohr 2004); and also subsamples of a complete data set can be used, for example by cross-checking data for X-ray or SZ proxies with weak lensing data (e.g., Vikhlinin et al. 2009). It has also been shown that characterising clusters using their gravitational potential instead of their mass can be an alternative way to reduce the cluster mass measurement bias (Tchernin et al. 2020).
A number of modified gravity (MG) theories (see, e.g., Koyama 2016) have been proposed over the past couple of decades in an attempt to better understand the as yet unexplained accelerated cosmic expansion. In these models, the law of gravity is typically modified, which means that clusters can become more (or less) massive than their counterparts in the standard Λ cold dark matter (ΛCDM) model. This can have consequences on the measured cluster number count, making the latter a potentially powerful cosmological probe to constrain deviations from GR.
A popular working example is the ( ) gravity model (see, e.g., Sotiriou & Faraoni 2010;De Felice & Tsujikawa 2010), which, like many other theories, predicts the presence of an extra scalar field that can lead to a strengthened force of gravity which is enhanced compared to Newtonian gravity. It is predicted that this could leave observational signatures in the formation of large-scale structures (LSS) that can make the model distinguishable from GR. A number of studies have therefore made use of LSS probes, such as the halo mass function probed by different proxies (see, e.g., Cataneo et al. 2015;Liu et al. 2016;Peirone et al. 2017), redshift-space distortions (e.g., Bose & Koyama 2017;Hernández-Aguayo et al. 2019;He et al. 2018), the cluster gas mass fraction (e.g., Li et al. 2016), the cluster temperature-mass relation (Hammami & Mota 2017) and SZ profile (De Martino 2016), and the clustering of clusters (Arnalte-Mur et al. 2017) to constrain these theories. For the current status and future prospect of testing gravity using galaxy clusters, see the recent reviews by Cataneo & Rapetti (2018); Baker et al. (2019).
The additional forces in the MG models can have complicated effects on the internal properties of galaxy clusters, including their density profile, dynamical mass and gas temperature. However, a robust and complete pipeline to account for these effects has not yet been implemented in tests that constrain MG theories using galaxy clusters. This paper is part of a series of works that aim to propose a framework for producing unbiased and self-consistent constraints of ( ) gravity (and other chameleon theories) using the halo mass function (and other LSS probes). The framework is described in Fig. 1 (for a more complete description, see Mitchell et al. 2018). Several parts of it are already complete. A simple yet powerful tanh fitting formula for the force enhancement was calibrated (Mitchell et al. 2018) using a suite of dark-matter-only simulations. This can be used to predict the enhancement of the dynamical mass of clusters. In particular, it can be used as an effective one-parameter description of the chameleon screening mechanism. This has already enabled us to calibrate a simple model of the ( ) enhancement of the halo concentration (Mitchell et al. 2019), and it is expected to simplify the calibration of the ( ) halo mass function in a future work. Our models for the enhancement of the dynamical mass and the concentration show a high level of accuracy for a wide and continuous range of background scalar field values, redshifts and halo masses. While these simplified and calibrated models of cluster properties were done using a particular ( ) model, their underlying logic implies that they will work for generic chameleon-type MG models, promising the generality of the framework.
In this work, we address the effects of ( ) gravity on cluster mass scaling relations. These are altered by the effects of the fifth force on the gravitational potential of haloes, which is intrinsically linked to the dynamical mass and gas temperature. According to , who employed a suite of non-radiative simulations to test the effects of ( ) gravity on a number of mass proxies, it is possible to map between scaling relations in ( ) gravity and GR using only the relation between the dynamical and true masses of a halo, which is accurately captured by our tanh fitting formula mentioned above (cf. Fig. 1). In this work, we test these predictions using nonradiative hydrodynamic simulations with much higher resolutions, and build on this by checking how the addition of full-physics effects such as cooling, star formation and feedback impact the accuracy. To do this, we make use of the first simulations that simultaneously incorporate both full-physics and ( ) gravity ). In addition, we propose and test a set of alternative mappings from the cluster mass scaling relations in GR to their ( ) counterparts, which again require only our tanh fitting formula.
The paper is organised as follows: in Sec. 2, several important aspects of the ( ) gravity theory are introduced, and we summarise the scaling relation mappings proposed by  and our novel alternative mappings; in Sec. 3, we describe our simulations and calculations of the halo masses and observable proxies; then, in Sec. 4, we present our results for the scaling relations of four mass proxies; finally, in Sec. 5, we summarise the results of this work and their significance for our framework.

BACKGROUND
We describe several key aspects of the Hu & Sawicki (2007) ( ) gravity model in Sec. 2.1. Then, in Sec. 2.2, we outline the effect of ( ) gravity on the dynamical mass of dark matter haloes and summarise predictions for how several cluster mass scaling relations are affected. We discuss the mappings between the ( ) and GR mass scaling relations proposed by   2) which is also tested in this work.

The ( ) gravity model
The ( ) gravity theory represents a modification to GR. The modification takes the form of the addition of a non-linear function, ( ), of the Ricci scalar curvature, , to the Einstein-Hilbert action : where is the determinant of the metric tensor, L M the Lagrangian density for matter and the (universal) gravitational constant. This alteration leads to the introduction of an extra tensor, , in the Einstein field equations: Here the tensors , and represent the Einstein tensor, the stress-energy tensor and the Ricci curvature, respectively. The d'Alembert operator is represented by , and ∇ is the covariant derivative associated with , the metric tensor. The scalar field ≡ d ( )/d represents the extra (scalar) degree of freedom of the theory. This mediates a so-called 'fifth force', which is an attractive force felt by massive particles. When this force is able to act, the effect is to enhance the strength of gravity by a factor of 4/3. The force can act on scales smaller than the Compton wavelength, C , which sets the physical range of it, where is the cosmological scale factor. Beyond C the amplitude of this force decays exponentially.  [Colour Online] Flow chart summarising the key steps of our framework to constrain modified gravity using galaxy clusters. For now, we focus on tests that use the halo mass function to constrain the ( ) gravity model proposed by Hu & Sawicki (2007). In previous works (Mitchell et al. 2018(Mitchell et al. , 2019, we used a suite of dark-matter-only simulations, which cover a wide range of mass resolutions and model parameters, to calibrate models for the enhancement of the dynamical mass (red dotted box) and the concentration (blue dotted box) of dark matter haloes in ( ) gravity. The model for the concentration can be used to convert between different halo mass definitions, which can allow for the derivation of the halo mass function in ( ) gravity using predictions from previous works (e.g., Cataneo et al. 2016). In this work (green dotted box), we use full-physics hydrodynamical simulations to test predictions that the dynamical mass enhancement can be used to convert scaling relations calibrated in ΛCDM into their ( ) gravity counterparts. These scaling relations can then be used to derive the cluster mass function, assuming ( ) gravity, starting from observational data. Finally, we will use Markov Chain Monte Carlo analysis to constrain the present-day background scalar field magnitude, | 0 |, using our theoretical predictions and observational calculations of the mass function (brown dotted box).
precision. For consistency with these tests, the chameleon screening mechanism is employed in the ( ) gravity model (e.g., Khoury & Weltman 2004a,b). It introduces an environment-dependent effective mass of the scalar field , which becomes more massive in high-density regimes. This has the effect of decreasing C and thus screening out the fifth force in high-density regions, including the solar system.
Following previous works, we focus on a representative variant of ( ) gravity proposed by Hu & Sawicki (2007, HS). This model is able to give rise to the late-time cosmic acceleration, while also showing consistency with solar system tests. It is characterised by the following prescription for ( ): where 2 ≡ 8¯M ,0 /3 = 2 0 Ω M , with 0 the Hubble constant, M,0 the present-day background matter density and Ω M the current dimensionless matter density parameter. The model has three parameters: 1 , 2 and . However, for a realistic set of cosmological parameters, it can be re-formulated to a good accuracy with just two free parameters: , which is fixed at 1 in this work, and the presentday background scalar field value 0 (we note that ( ) and 0 represent background values for the remainder of this paper). We investigate models with values | 0 | = 10 −6 and | 0 | = 10 −5 , referring to these as F6 and F5, respectively. F5 represents a stronger modification to GR than F6, and allows higher-mass haloes to be unscreened.

Mass scaling relations in ( ) gravity
Various techniques can be employed to measure the mass of clusters. One approach is to use the gravitational lensing of background sources. In ( ) gravity, apart from an overall rescaling of the lens mass by (1 + ) −1 , which is very close to unity for realistic choices of the ( ) parameters, photons are unaffected by the fifth force, meaning that the lensing effect of a given cluster is the same as in GR. In this work, we refer to the mass inferred from lensing as the 'true' mass, true , and this satisfies the following: The mass can alternatively be measured using probes such as the X-ray temperature. These can provide a measure of the 'dynamical' mass, dyn , which we define as the mass that is felt by nearby massive test particles. Because the gravitational potential of unscreened haloes is enhanced in ( ) gravity, the dynamical mass of haloes is also enhanced by the same factor: true . On the other hand, in GR, the dynamical mass is expected to be equal to the true mass: GR dyn = GR true = GR . Using a suite of -body simulations covering a wide range of resolutions, Mitchell et al. (2018) showed that the enhancement of the dynamical mass in ( ) gravity can be accurately described by a simple tanh formula, given by, The parameter 1 was found to be approximately constant with value 2.21 ± 0.01, while a simple power law, whose slope was motivated by theory (which predicts it to be 3/2), was calibrated for 2 : 2 = (1.503 ± 0.006) log 10 | ( )| 1 + + (21.64 ± 0.03).
This depends on just a single parameter, | |/(1+ ), and provides a powerful yet simple approach to describe the chameleon screening mechanism in ( ) gravity.
In observational studies, obtaining detailed and high-quality X-ray and spectral data for determination of the cluster dynamical mass is often an expensive task, requiring long exposure times. It is therefore more common to predict the cluster mass using its relationship with more easily measured observables, including the temperature, the X-ray luminosity, and the Compton -parameter of the SZ effect, and its X-ray analogue. These mass 'proxies' have a one-to-one mapping with the mass because of the link between the gravitational potential of a cluster and its temperature. During cluster formation, baryonic matter is accreted onto the dark matter halo from its surroundings. The gravitational potential energy of the gas is converted into kinetic energy as it falls in. During accretion, the in-falling gas undergoes shock heating, resulting in the conversion of its kinetic energy into thermal energy. The resulting self-similar model for cluster mass scaling relations predicts that the gravitational potential alone can determine the thermodynamical properties of a cluster (Kaiser 1986).
In this work, we study the effects of HS ( ) gravity on the scaling relations for four mass proxies: the gas temperature gas , the X-ray luminosity X , the integrated SZ flux, given by the Compton -parameter SZ , and its X-ray analogue X . The X-ray luminosity within radius from the cluster centre is given by, where gas ( ) is the gas density profile. The SZ parameter is related to the integrated electron pressure of the cluster gas, and is given by, where T is the Thomson scattering cross section, e is the electron rest mass, is the speed of light and e is the electron number density. Meanwhile, the X parameter (Kravtsov et al. 2006) is equivalent to the product of the gas mass and the mass-weighted gas temperature,¯g as (< ): It has been shown in previous studies (e.g., Fabjan et al. 2011) that X and SZ have comparatively low scatter as mass proxies and are relatively insensitive to dynamical processes including cluster mergers. It has also been found that their scaling relations with the mass show good agreement with the self-similar model predictions even after the inclusion of full-physics effects such as feedbacks, which can heat up and blow gas out from the central regions (e.g., Fabjan et al. 2011;Truong et al. 2018;Cui et al. 2018).
Because these cluster mass proxies are closely related to the dynamical mass of clusters, and because in ( ) gravity the dynamical mass can be cleanly modelled (see above), it is natural to expect that cluster observable-mass scaling relations in ( ) gravity can be modelled given their counterparts in GR. To study the effects of ( ) gravity on the scaling relations for these proxies, we adopt two methods which are described in the sections that follow.

Effective density approach
The ( ) gravity effective density field, eff , was originally defined by He et al. (2014): where is the true density field, corresponding to the intrinsic mass of simulation particles. Written in terms of the effective density, the Poisson equation of ( ) gravity is cast into the same form as in Newtonian gravity: where is the total gravitational potential, including contributions from the fifth force. It follows that the mass of haloes computed using the effective density field is equivalent to the dynamical mass.
Using non-radiative simulations run for the F5 model and GR, He & Li (2016) generated halo catalogues using the effective density field. The radius, eff 500 , of these haloes enclosed an average effective density of 500 times the (true) critical density of the Universe. Both the total true and dynamical mass were computed within this radius, and the cluster observables were computed using all enclosed gas particles.
Analysing these data,  found that the mapping between the gas temperature and the dynamical mass is identical for haloes in GR and ( ) gravity: This result is in line with the self-similar model, since two haloes with the same dynamical mass (which we recall has also been computed within the same radius), GR = ( ) dyn = dyn , also have the same gravitational potential, = ( dyn )/ eff 500 . The authors also found that, outside the core region, the gas density profiles of haloes in GR are enhanced by a factor ( ) dyn / ( ) true with respect to haloes in ( ) gravity which have the same dynamical mass. This is because the gas density profile follows the true density profile more closely than the effective density profile, which itself is a result of the fact that clusters form from very large regions in the Lagrangian space, so that the ratio between the baryonic and total masses within clusters resembles the cosmic mean, Ω B /Ω M (White et al. 1993), in which Ω B is the present-day baryonic density parameter. The extra forces in MG theories and feedbacks from galaxy formation can add further complications to this through their effects on the gas density profiles, especially in the inner regions; however, as we will show in the following paragraph, the good agreement between the GR and rescaled ( ) gas density profiles still holds in the outer halo regions.
We have replicated the procedure adopted by He & Li (2016) using our full-physics and non-radiative simulations (for full details, see Sec. 3). In Fig. 2, the stacked temperature and gas density profiles of haloes from mass bins 10 13.7 < dyn (< eff 500 ) < 10 14.0 and 10 13.0 < dyn (< eff 500 ) < 10 13.3 are shown in the first and second columns from the left, respectively. The radial range is shown up to the mean logarithm of eff 500 (which is almost exactly the same for GR, F6 and F5). For the non-radiative temperature profiles, shown in the third row, it is clear that the F6 and F5 predictions agree very well with GR in the outer regions. There is also encouraging agreement for the full-physics data, although there is a small disparity between F5 and GR in the outer regions for the higher-mass bin. For the ( ) gravity gas density profiles, the results both with and without the  full-physics true Figure 2. [Colour Online] Median gas density profiles (top two rows) and median temperature profiles (bottom two rows) of FOF groups from two mass bins: 13.7 < log 10 ( / ) < 14.0 (high-mass) and 13.0 < log 10 ( / ) < 13.3 (low-mass). The group data from the non-radiative and full-physics SHYBONE simulations (see Sec. 3.1) has been used. In addition to GR (red solid lines), the profiles for F6 (blue lines) and F5 (green lines) are shown. Rescaled ( ) gravity profiles (dashed lines) are shown, along with the unaltered profiles (dotted curves). The rescalings correspond to the effective density (left two columns) and true density (right two columns) approaches discussed in Sec. 2.2. For the effective approach, the maximum radius shown is eff 500 (see Sec. 2.2.1) and the halo mass is the total dynamical mass within this radius. For the true approach, the maximum radius shown is true 500 (see Sec. 2.2.2) and the halo mass is the total true mass within this radius.
shown. As was found by , the rescaled ( ) gravity profiles (shown by the dashed curves) agree very well with GR in the outer regions. Again, there is also promising agreement for the full-physics profiles.
These results suggest that for haloes in ( ) gravity and GR which have the same dynamical mass, GR = ( ) dyn , the follow-ing relation is expected to apply: where and represent indices of power, and we note that gas represents the intrinsic (not effective) gas density. Using Eqs. (8)-(10) for the mass proxies, this relation can be applied to derive the following mappings between the respective mass scaling relations in GR and ( ) gravity: Note that to obtain these relations, the two integrations in Eq. (14) have used the same upper limit, = eff 500 , for GR and ( ) gravity, as mentioned above.  demonstrated an accuracy ≈ 3% for the SZ and X relations and ≈ 13% for X . These quantities are all cumulative: they are computed by summing over the entire volume of a halo up to some maximum radius (in this case eff 500 ). Cumulative quantities typically depend more on the outer regions, which contain most of the volume and mass, than on the inner regions. Therefore, even though the ( ) gravity profiles (with appropriate rescaling) do not agree with GR for the inner regions (see Fig. 2), this is expected to have a negligible contribution overall to these mass proxies and explains why  found such a high accuracy for these mappings.
In our work, we will expand on these tests by using full-physics hydrodynamical simulations to check how the addition of effects such as cooling and feedback can alter the accuracy of the mappings defined by Eqs. (15)-(17) and the temperature equivalence given by Eq. (13). Our tests with the non-radiative runs can also be used as a check for consistency with , who used a different simulation code and ( ) gravity solver.

True density approach
Mappings can also be done for haloes identified with the true density field. For these haloes, the radius, true 500 , would enclose an average true density of 500 times the critical density of the Universe. For haloes in ( ) gravity and GR with the same true mass, GR = ( ) true , the total gravitational potential at true 500 in ( ) gravity would be a factor of ( ) dyn / ( ) true higher than in GR (where both the dynamical and true mass are measured within true 500 ). According to the self-similar model predictions, the gas temperature in ( ) gravity is expected to be higher by the same factor: On the other hand, for haloes in ( ) gravity and GR with the same true mass, the gas density profiles are expected to agree in the outer regions.
To check these assumptions, let us once again examine Fig. 2, this time looking at the third and fourth columns from the left, which show the stacked gas density and temperature profiles for groups in mass bins 10 13.7 < true < true 500 < 10 14.0 and 10 13.0 < true < true 500 < 10 13.3 , respectively. The radial range is shown up to the mean logarithm of true 500 for all profiles. Referring to the bottom two rows, which show the non-radiative and full-physics temperature profiles, it appears that the ( ) gravity profiles with the ( ) dyn / ( ) true rescaling (shown by the dashed curves) show reasonable agreement with GR in the outer regions. Again, the only exception is for the high-mass bin of the full-physics data, where there is a small disparity between F5 and GR. Looking at the top two rows, the ( ) gravity gas density profiles agree very well with GR in the outer-most regions for both mass bins.
These results for the temperature and gas density profiles yield the following predictions for haloes in GR and ( ) gravity with GR = ( ) where this time the two integrations both have upper limit = true 500 . This prediction yields the following new mappings between the mass scaling relations in ( ) gravity and GR: For the SZ and X mappings, the ( ) dyn / ( ) true factor comes from the dependence on the gas temperature to power one in Eqs. (9) and (10). On the other hand, the X-ray luminosity, given by Eq. (8), depends on the gas temperature to power half, which means that the corresponding ( ) gravity and GR scaling relations are expected to differ by factor only. In Sec. 4, we show the results of our tests of these alternative predictions using both our non-radiative and full-physics simulations.
In this section, we have referred to two different definitions of the halo radius: eff 500 and true 500 . The radius eff 500 is defined in terms of the effective density field. For an unscreened halo in ( ) gravity, the effective density is up to 4/3 times greater than the true density. As such, eff 500 is typically a higher radius than true 500 for haloes in ( ) gravity. On the other hand, the two radii are exactly the same in GR, where the effective density field is equivalent to the true one.

SIMULATIONS AND METHODS
In Sec. 3.1, we describe the non-radiative and full-physics simulations that are used in this work. Then, in Sec. 3.2, we describe how we have measured the cluster mass and four observable mass proxies from these simulations.

Simulations
The results discussed in this work were generated using a subset of the SHYBONE simulations ). These full-physics hydrodynamical simulations employ the IllustrisTNG galaxy formation model (Weinberger et al. 2017;Pillepich et al. 2018) and include runs for both GR and HS ( ) gravity. The simulations have been run using the code (Springel 2010), which is a highly parallel and optimised code for hydrodynamical cosmological simulations. The code features a modified gravity solver which uses adaptive mesh refinement to accurately measure the fifth force in high-density environments. For every full-physics run used in this work, we also utilise a non-radiative counterpart which does not include cooling, star formation or stellar and black hole feedback processes.
Both the full-physics and non-radiative simulations span a comoving box of length 62ℎ −1 Mpc. These runs each start with 512 3 dark matter particles and the same number of initial gas resolution elements ( and an average gas cell mass of gas ≈ 2.5 × 10 7 ℎ −1 . In addition to GR, the runs include the F6 and F5 HS models, all starting from identical initial conditions at = 127. In the calculation of the gas temperature, we have assumed that the primordial hydrogen mass fraction has a value H = 0.76 and set the adiabatic index to = 5/3 (for a monatomic gas). For the non-radiative simulations we assume that the gas is made up of fully ionised hydrogen and helium.

Group catalogues
The group catalogues were constructed using the code implemented in AREPO, which employs a standard friends-of-friends (FOF) algorithm combined with an un-binding method to identify the bound structures within a FOF group and the gravitational potential minimum of the objects (Springel et al. 2001) using the true density field. For each group, both radii true 500 and eff 500 were computed around the gravitational potential minimum, enclosing, respectively, average true and effective densities of 500 times the critical density of the Universe. For each radius definition, the total enclosed dynamical and true masses were measured, in addition to the group observables. Quantities measured within eff 500 have been used to test the scaling relation mappings from the effective density approach described in Sec. 2.2.1, while the quantities measured within true 500 have been used to test the predictions of the true density approach discussed in Sec. 2.2.2.
In the computation of the halo temperature, we have excluded the core regions in which the complex thermal and dynamical processes during cluster formation and evolution can lead to a significant degree of dispersion between the halo temperature profiles. We have set the core region to the radial range < 0.15 , where can be either eff 500 or true 500 . This range is consistent with previous studies of cluster scaling relations (e.g., Fabjan et al. 2011;Truong et al. 2018).
The halo gas temperature has been computed using a mass-weighted average: where gas, and are, respectively, the mass and temperature of gas cell . The summations have been performed over all gas cells whose positions fall within the radial range 0.15 < < . The integrated SZ flux is given by: where e, is the number of electrons in gas cell and the sum includes the same cells as for¯g as . The X-ray analogue of the integrated SZ flux is equal to the product of the total gas mass gas , of all gas cells within 500 , and¯g as : Finally, the X-ray luminosity is calculated using: where gas, is the gas density of gas cell and the summation is performed over the same gas cells as for the¯g as calculation.

RESULTS
In Sec. 4.1, we discuss our results for the cluster scaling relations in HS ( ) gravity. Then, in Sec. 4.2, we test the validity of our analytical tanh formula for the dynamical mass enhancement, given by Eq. (6), in the presence of full physics. There we will also present an example in which we map between the GR and ( ) mass scaling relations based on this approximate fitting formula, rather than the actual values of ( ) dyn / ( ) true from the simulations.

Scaling relations
Using our simulation data, we have tested the scaling relation mappings described in Sec. 2.2. Due to the small box-size 62ℎ −1 Mpc (comoving) of our simulations, there are only a few haloes with mass greater than 10 14 . We therefore show only the range 10 13 to 10 14 in Figs. 3-6. This range typically includes ∼ 100 haloes for a given model. For F5, the haloes are typically unscreened throughout this range, while for F6 the haloes are either partially screened (at lower mass) or completely screened (at higher mass). In addition to showing data points for individual haloes, all plots include curves showing a moving average. This is computed using a moving window of fixed size equal to 10 haloes, for which the mean logarithm of the mass and the median proxy are displayed. Plots of the smoothed relative difference between the ( ) gravity curves and the GR curves are also displayed.
For the panels labelled 'effective' in Figs. 3-6, all measurements of the mass and observables have been taken within eff 500 (see Sec. 2.2.1), and the relations are plotted against the dynamical mass. This allows the effective density mappings given by Eqs. (13) and (15)-(17), originally proposed by , to be tested. On the other hand, an outer radius true 500 (see Sec. 2.2.2) is imposed for all measurements for the data displayed in the panels labelled 'true'. These are plotted against the true mass, and can be used to test the true density mappings given by Eqs. (18) and (20)   Gas temperature plotted as a function of the mass for FOF groups from the non-radiative and full-physics SHYBONE simulations (see Sec. 3.1). The curves correspond to the median mass-weighted temperature and the mean logarithm of the mass computed within a moving window of fixed size equal to 10 haloes. Data has been included for GR (red solid lines) together with the F6 (blue lines) and F5 (green lines) ( ) gravity models. Rescalings to the ( ) gravity temperature have been carried out as described in Secs. 2.2.1 and 2.2.2. For the 'true' density approach, the rescaled data (dashed lines) is shown along with the unaltered data (dotted lines). For this data, the mass corresponds to the total true mass within the radius true 500 , and the temperature has also been computed within this radius. For the 'effective' density approach, no rescaling is necessary, the mass corresponds to the total dynamical mass within eff 500 , and the temperature has also been computed within this radius. Data points are displayed, with each point corresponding to a GR group (red points) or to a group in F6 (blue points) or F5 (green points), including the rescaling for the 'true' density data. Bottom row: the smoothed relative difference between the ( ) gravity and GR curves in the above plots.

Temperature scaling relations
The results for the gas temperature scaling relations are shown in Fig. 3. The non-radiative data is displayed in the left two columns and the full-physics data is shown in the right two columns. For all models and hydrodynamical schemes, the data follows a power-law behaviour, as expected from the self-similar model. The correlation is particularly tight for the non-radiative runs. These contain gas and dark matter particles, but do not feature baryonic processes such as cooling, stellar and black hole feedback and star formation. It is therefore expected that the thermodynamical properties can be largely determined from the gravitational potential, which is observed in the results. On the other hand, there is a greater amount of scatter in the halo data for the full-physics runs, and the gas temperature is typically higher. This can be explained by the inclusion of feedback mechanisms which act as a further source of heating of the surrounding gas and cause some departures from self-similarity. For the effective density approach, Eq. (13) is expected to hold: the temperature is predicted to be equal for haloes in GR and ( ) gravity with the same dynamical mass. In Fig. 3, the non-radiative and full-physics results from our effective catalogue are shown in the first and third columns from the left, respectively. For both F6 and F5, there is excellent agreement with the GR data, with typical agreement 5%. This agreement for the non-radiative data backs up the findings from , while the full-physics results do not show clear evidence for a departure from Eq. (13) caused by feedback processes and cooling. These results agree with the selfsimilar model predictions: two haloes in ( ) gravity and GR which have the same dynamical mass dyn (and therefore the same radius eff 500 ) also have the same gravitational potential, = dyn / eff 500 . In order to test the new mappings predicted by the true density approach, the temperature of each halo in ( ) gravity has been divided by the mass ratio  Compton -parameter of the SZ effect, plotted as a function of the mass for FOF groups from the non-radiative and full-physics SHYBONE simulations (see Sec. 3.1). The curves correspond to the median SZ versus the mean logarithm of the mass computed within a moving window of fixed size equal to 10 haloes. Data has been included for GR (red solid lines) together with the F6 (blue lines) and F5 (green lines) ( ) gravity models. Rescalings to SZ have been carried out as described in Sec. 2.2, and both the rescaled (dashed lines) and unaltered (dotted lines) results are shown. For the 'effective' density approach (see Sec. 2.2.1), the mass corresponds to the total dynamical mass within eff 500 , and SZ has also been computed within this radius. For the 'true' density approach (see Sec. 2.2.2), the mass corresponds to the total true mass within the radius true 500 , and SZ has also been computed within this radius. Data points are displayed, with each point corresponding to a GR group (red points), or to a group in F6 (blue points) or F5 (green points) with the relevant rescaling applied to SZ . Bottom row: the smoothed relative difference between the ( ) gravity and GR curves in the above plots. the non-radiative results in the second column, this indeed appears to be the case, with an excellent agreement that is generally within just a few percent. For the full-physics data there is still reasonable 10% agreement, but the F5 temperature appears to be lower than the GR temperature for log 10 true < true 500 / 13.5. This small deviation is consistent with the full-physics temperature profiles shown in Fig. 2. The plots in the bottom right of that figure show the temperature profiles with the ( ) dyn / ( ) true rescaling applied. The profiles are shown for the mass bins 10 13.7 < true < true 500 < 10 14.0 and 10 13.0 < true < true 500 < 10 13.3 , which correspond to the high-and low-mass end of the mass range shown in Fig. 3. For the high-mass full-physics profiles, the rescaled F5 profile is clearly lower than the profile in GR across most of the radial range. This can explain the lower rescaled F5 temperature observed at the high-mass end of the full-physics data. On the other hand, for the lower-mass bin with full-physics and for the non-radiative data the agreement between the rescaled ( ) gravity and GR temperature profiles is very good, particularly at the outer radii which have greater overall contribution to the mass-weighted temperature. Similar agreement is shown between the ( ) gravity and GR profiles for the plots in the bottom-left of Fig. 2. Again, for the high-mass full-physics data there is some deviation between F5 and GR, particularly in the outermost regions. But this is not as noticeable as for the profiles with the true density rescalings.
The small difference between F5 and GR is likely to be caused by a difference in the levels of feedback -which itself is determined by the interrelations between modified gravity (including screening or the lack of it) and baryonic physics -in these higher-mass haloes for the two models. It is encouraging that this results in just a small effect overall.

SZ and X scaling relations
Our results for the SZ and X scaling relations are shown in Figs. 4 and 5. The SZ and X parameters are, by definition, tightly correlated. Their results therefore follow similar patterns. They show very [Colour Online] X-ray analogue of the Compton -parameter plotted as a function of the mass for FOF groups from the non-radiative and full-physics SHYBONE simulations (see Sec. 3.1). The curves correspond to the median X versus the mean logarithm of the mass computed within a moving window of fixed size equal to 10 haloes. Data has been included for GR (red solid lines) together with the F6 (blue lines) and F5 (green lines) ( ) gravity models. Rescalings to X have been carried out as described in Sec. 2.2, and both the rescaled (dashed lines) and unaltered (dotted lines) results are shown. For the 'effective' density approach (see Sec. 2.2.1), the mass corresponds to the total dynamical mass within eff 500 , and X has also been computed within this radius. For the 'true' density approach (see Sec. 2.2.2), the mass corresponds to the total true mass within the radius true 500 , and X has also been computed within this radius. Data points are displayed, with each point corresponding to a GR group (red points), or to a group in F6 (blue points) or F5 (green points) with the relevant rescaling applied to X . Bottom row: the smoothed relative difference between the ( ) gravity and GR curves in the above plots.
tight correlations with the halo mass, with the halo data following a narrower band than the temperature and X-ray luminosity. There is more scatter in the full-physics data, but there are no clear outliers, unlike for the temperature and the X-ray luminosity data. This is because of the competing effects of feedback processes on the gas density and gas temperature (Fabjan et al. 2011). Comparing the non-radiative and full-physics profiles in Fig. 2, it can be seen that the additional processes in the full-physics runs cause haloes to have a lower gas density, particularly at the inner regions, and a higher gas temperature. This is caused by stellar and black hole feedbacks, which generate high-energy winds that heat up the surrounding gas and blow it out from the central regions. Such competing effects are approximately cancelled out in the product of the gas density with the gas temperature, as in Eqs. (9) and (10).
In order to test the mappings predicted by the effective density approach, given by Eqs. (15) and (16), the SZ and X values measured for ( ) gravity have been multiplied by the mass ratio ( ) dyn / ( ) true . For the non-radiative plots in Figs. 4 and 5, there is excellent agreement between this rescaled data and GR. There is also a strong agreement for the higher-mass full-physics data. For the mass range 13.2 < log 10 ( dyn (< eff 500 ) −1 ) < 13.5, however, there is some disparity of 20% between the rescaled F5 data and GR.
A similar level of accuracy is observed for the mappings given by Eqs. (20) and (21), which are predicted by the true density approach. To test these, SZ and X are divided by ( ) dyn / ( ) true to generate rescaled data for F6 and F5. Again, this data shows excellent agreement, within a few percent, with GR for the non-radiative simulations and the high-mass end of the full-physics data. But for the lower-mass full-physics data there is a significant disagreement between F5 and GR of up to ∼ 30%, which is higher than for the effective density rescalings.
The disparities found in the low-mass full-physics data can be explained using Fig. 2. Looking at the full-physics profiles for the true density mass bins, it is observed that for a large portion of the inner halo regions the gas density is higher in F5 than in GR. These profiles converge at ≈ 10 2.5 kpc for both mass bins. The high-mass haloes have higher overall radius true 500 , which means [Colour Online] X-ray luminosity plotted as a function of the mass for FOF groups from the non-radiative and full-physics SHYBONE simulations (see Sec. 3.1). The curves correspond to the median luminosity versus the mean logarithm of the mass computed within a moving window of fixed size equal to 10 haloes. Data has been included for GR (red solid lines) together with the F6 (blue lines) and F5 (green lines) ( ) gravity models. Rescalings to the luminosity have been carried out as described in Sec. 2.2, and both the rescaled (dashed lines) and unaltered (dotted lines) results are shown. For the 'effective' density approach (see Sec. 2.2.1), the mass corresponds to the total dynamical mass within eff 500 , and the luminosity has also been computed within this radius. For the 'true' density approach (see Sec. 2.2.2), the mass corresponds to the total true mass within the radius true 500 , and the luminosity has also been computed within this radius. Data points are displayed, with each point corresponding to a GR group (red points), or to a group in F6 (blue points) or F5 (green points) with the relevant rescaling applied to the luminosity. Bottom row: the smoothed relative difference between the ( ) gravity and GR curves in the above plots. that the profiles are converged for a large portion of the outer radii. This means the disparities at lower radii have a negligible overall contribution to the integrals for SZ and X . But for the lower-mass haloes, the profiles are converged only for a small portion of the overall radius range, causing SZ and X to be greater in F5 than in GR. A similar reasoning can be used to explain the disparities for the results with the effective density rescaling, although the difference in agreement at the outer regions for each mass bin is not quite as substantial here, which explains why the full-physics data for the effective density approach shows less overall deviation between F5 and GR in Figs. 4 and 5.
The difference between the full-physics ( ) and GR scaling relations at low mass is likely to be explained by baryonic processes such as feedback which are absent in the non-radiative simulations. However, in studies of clusters, these lower-mass groups are of less interest. The strong agreement at the higher masses is therefore very encouraging for our framework to constrain ( ) gravity using the high-mass end of the halo mass function.

X-ray luminosity scaling relations
Our results for the X-ray luminosity are shown in Fig. 6. Compared with the temperature, SZ and X data, the X-ray luminosity is much more scattered, with large dispersion in lower-mass haloes in particular. One explanation for this is that the X-ray luminosity, defined in Eq. (8), depends on the gas density to power two. This means that the inner regions of the group, which have a higher gas density than the outer regions, have a greater overall contribution to the X integral than for the other observables discussed in this work. The inner halo regions are expected to be more impacted by unpredictable dynamical processes during cluster formation, including halo mergers. In particular, they are more prone to gas heating and blowing-out of gas caused by feedback mechanisms. While the competing effects of these processes on the gas density and gas temperature roughly cancel for the SZ and X observables, this is not the case for X , which depends on the gas density to power two and the gas temperature to power half. This results in a number of significant outliers, as can be seen in the full-physics data of Fig. 6.
For the mapping defined using the effective density field, given by Eq. (17), it has been expected that the GR X-ray luminosity should be equal to the ( ) gravity value multiplied by the factor ( ) dyn / ( ) true 2 . From Fig. 6, the rescaled data in F5 appears to be higher than in GR by ∼ 30% on average for both the non-radiative and the full-physics simulations. A similar level of deviation is also observed for the true density results, where Eq. (22) predicts that the GR and ( ) gravity X-ray luminosity should be equal after the values in ( ) gravity are divided by the factor ( ) dyn / ( ) true 1/2 . Again, the rescaled F5 X-ray luminosity is significantly greater than in GR on average.
As for the SZ and X mappings, the disparity found here can be explained by looking at the gas density profiles in Fig. 2. For both the non-radiative and full-physics data, the gas density in the inner halo regions is greater for F5 (with appropriate rescaling applied) than for GR. Because the inner regions have a greater contribution to the Xray luminosity than for other proxies, as described above, this causes these differences in the inner regions to become significant overall, even for the non-radiative data for which the F5 and GR profiles are converged above a lower radius. This results in the general offset for the full range of masses as shown in Fig. 6. As described above, the X-ray luminosity is also more strongly influenced by feedback processes, which can further increase the offset between F5 and GR if the feedback behaves differently in these two models.
Our observation that the mappings have a poorer performance for the X-ray luminosity than for the other proxies is consistent with , who observed a disparity of ∼ 13%. Unless further corrections are applied to account for the unpredictable effects of feedback in the mappings, therefore, the X-ray luminosity is unlikely to be a reliable proxy for mass determination in accurate cluster tests of ( ) gravity.

Halo mass ratio calculation
For the results discussed in Sec. 4.1, the rescalings to the observables in ( ) gravity have been computed using direct measurements of the true mass and the dynamical mass from the simulations. However, for studies of clusters using real observations, measurements of both the true mass and dynamical mass are unlikely to be available. In this case, our analytical model for the ratio of the dynamical mass to the true mass, given by Eq. (6), can be used. This formula was calibrated by Mitchell et al. (2018) using a suite of dark-matter-only simulations and tested for wide and continuous ranges of scalar field values, redshifts and halo masses.
To check how this model performs for data that includes full physics, we have plotted this on top of actual measurements of the dynamical mass enhancement for the FOF groups in the SHYBONE simulations. This is shown in Fig. 7. Data for both the effective density and true density catalogues have been included, for which the dynamical and true halo masses have been measured within eff 500 and true 500 , respectively. We have made use of all available data with true > 10 13 , including haloes with true ∼ 10 14.5 . These results indicate that there is very good agreement between the analytical predictions and the actual data for both F6 and F5, regardless of the hydrodynamical scheme that is employed. Interestingly, even though Eq. (6) was originally calibrated using measurements of the dynamical and true mass within true 500 , it still performs very well for data measured within eff 500 , which is typically a higher radius. We have also tested the mappings of Eqs. (16) and (21) for the X parameter, with Eq. (6) used to compute the required rescalings to the ( ) gravity data. This is shown in Fig. 8. From comparing this plot with Fig. 5, it can be seen that there is almost no difference in the rescaled data in both figures. This confirms that Eq. (6) can be applied to derive the mappings between GR and ( ) scaling relations for, at least, the mass range 10 13 < 500 < 10 14 . Given the very good agreement up to 10 14.5 shown in Fig. 7, it is expected that our formula can be applied at higher masses as well.

SUMMARY, DISCUSSION AND CONCLUSIONS
Accurate determination of the cluster mass is a vital requirement of tests that make use of LSS probes including cluster counts to constrain cosmological models. It is particularly important in tests of chameleon-type MG theories, in which the additional forces present can alter the internal properties of clusters, including the concentration, the dynamical mass and the temperature. This means that scaling relations between the cluster mass and observable proxies that have been determined assuming GR are unlikely to apply in a universe with, for example, chameleon ( ) gravity.
This work is part of a series (Mitchell et al. 2018(Mitchell et al. , 2019 which aims to create the first complete and robust pipeline for accounting for the effects of ( ) gravity on the internal properties of galaxy clusters. Highlights so far include the calibration of analytical models that can accurately describe the enhancement of the dynamical mass (Mitchell et al. 2018) and concentration (Mitchell et al. 2019) of haloes in HS ( ) gravity. These models are accurate for a wide and continuous range of masses, redshifts and values of the background scalar field. In particular, we have found a powerful way to describe the chameleon screening mechanism with just a single parameter: | |/(1 + ), and in a future work will use this to calibrate a general model for the halo mass function in ( ) gravity.
In this work, we have made use of the first full-physics simulations that have been run for both GR and ( ) gravity (along with non-radiative counterparts), to study the effects of the fifth force of ( ) gravity on the scaling relations between the cluster mass and four observable proxies: the gas temperature (Fig. 3), the SZ and X parameters (Figs. 4 and 5) and the X-ray luminosity (Fig. 6). To understand these effects in greater detail, we have also examined the effects of both ( ) gravity and full-physics on the gas density and temperature profiles (see Fig. 2). In doing so, we have been able to test two methods for mapping between scaling relations in ( ) gravity and GR.
The first method was proposed by . This proposes a set of mappings, given by Eqs. (13) and (15)-(17), that can be applied to haloes whose mass and radius are measured using the effective density field (see Sec. 2.2.1). A second, new, approach is proposed in Sec. 2.2.2, and predicts another set of mappings, given by Eqs. (18) and (20)-(22), that can be applied to haloes whose mass and radius are measured using the true density field. Both sets of mappings involve simple rescalings that depend only on the ratio of the dynamical mass to the true mass in ( ) gravity. As shown by Figs. 7 and 8, even with the inclusion of full-physics processes this ratio can be computed with high accuracy using our analytical tanh formula, which is given by Eq. (6).
For the mass-weighted gas temperature and the SZ and X observables, we found that the F6 and F5 scaling relations, with appropriate rescaling applied (using either method discussed above), match the GR relations to within a few percent for the full massrange tested for the non-radiative simulations. With the inclusion of full-physics effects such as feedbacks, star formation and cooling, the rescaled SZ and X scaling relations continue to show excellent The data points correspond to FOF groups from the non-radiative and full-physics SHYBONE simulations (see Sec. 3.1 for details). Data is included for F6 (blue) and F5 (green), which are the HS ( ) gravity models with = 1 and present-day background scalar field | 0 | = 10 −6 and | 0 | = 10 −5 , respectively. For the data labelled 'effective', the dynamical and true mass have been measured within the radius eff 500 (defined in Sec. 2.2.1), while the radius true 500 (defined in Sec. 2.2.2) has been used for the data labelled 'true'. Analytical predictions for the mass enhancement have been computed using Eq. (6) and are shown (dashed curves) for each model.  Fig. 5, however the rescalings of X in the ( ) gravity data have been generated using an analytical tanh formula, given by Eq. (6). agreement with GR for mass 500 10 13.5 , which includes group-and cluster-sized objects. These proxies also show relatively low scatter as a function of the cluster mass, compared with other observables. SZ and X are therefore likely to be suitable for accurate determination of the cluster mass in tests of ( ) gravity. The mappings for the gas temperature show a very high accuracy for lower-mass objects, but show a small 5% offset between F5 and GR for higher-mass objects.
The mappings do not work as well for the X-ray luminosity X , for which the F5 relations after rescaling are typically enhanced by ∼ 30% compared with GR. This is caused by the unique dependency of X on the gas density to power two, and the gas temperature to power half, which means that the inner halo regions have a greater contribution than for the other proxies and the competing effects of feedback on the temperature and gas density profiles are less likely to cancel out. This issue, in addition to the fact that X has a highly scattered correlation with the cluster mass, means that this proxy is unlikely to be suitable for cluster mass determination in tests of ( ) gravity. We note that the box size 62ℎ −1 Mpc of the simulations used in this work is more suited to studying galaxy-sized objects than group-or cluster-sized objects. Indeed, there are only ∼ 100 objects with 500 > 10 13 and ∼ 5-10 objects with 500 > 10 14 in the simulations. This makes it impossible to test the mappings discussed in this work for the most massive galaxy clusters to be observed. We are therefore currently preparing to run larger simulations, which employ a re-calibrated full-physics model, that can be used to reliably probe halo masses up to 500 ∼ 10 15 . We will leave the analysis of these simulations to a future work.