Testing Screening Mechanisms with Mass Profiles of Galaxy Clusters

We present \textsc{MG-MAMPOSSt}, a license-free code to constrain modified gravity models by reconstructing the mass profile of galaxy clusters with the kinematics of the cluster's member galaxies. We describe the main features of the code and we show the capability of the method when the kinematic information is combined with lensing data. We discuss recent results and forecasts on two classes of models currently implemented in the code, characterized by different screening mechanisms, namely, chameleon and Vainshtein screening. We further explore the impact of possible systematics in view of application to the data from upcoming surveys. This proceedings summarizes the results presented at the ALTECOSMOFUN workshop in September 2021.


Introduction
Scalar-tensor theories of gravity represent one of the most general classes of modified gravity (MG) models viable at cosmological scales, characterized by a quite large number of interesting features (see, e.g., [1]). This family of theories includes several subclasses, such as the Horndeski (e.g., [2]) and the more general Degenerate Higher-Order Scalar-Tensor (DHOST hereafter) theories (e.g., [3,4]).
A new dynamical scalar field propagates an additional fifth force which enhances gravity locally on sufficiently large scales, producing possible detectable imprints on the formation and evolution of cosmological structures. In order to match observational constraints on small scales and dense environments, the fifth force should be further suppressed through a so-called screening mechanism, which relies on the non-linear interaction of the scalar field (e.g., [5]).
In this framework, galaxy clusters constitute a powerful tool to investigate modification of gravity at large scales (see [6] for a recent review); in particular, with the combination of the cluster's mass profiles derived with lensing and internal kinematic analyses of the hot intra-cluster gas (e.g., [7][8][9]) and of cluster's member galaxies (e.g., [10,11]), it is possible to constrain departures from General Relativity (GR) in a complementary way with respect to other cosmological and astrophysical observations.
To this end, we developed MG-MAMPOSST, a FORTRAN code which determines cluster mass profiles in modified gravity from the kinematics of cluster members; the program is based upon the MAMPOSST method of [12]. While the latter assumes that the gravitational interaction is described by GR, MG-MAMPOSST explores two popular and quite general classes of MG models, characterized by distinct screening mechanisms: chameleon screening (CS) and Vainshtein screening (VS).
The code further implements a lensing simulation to investigate the capability of a joint lensing and kinematic analysis to constrain the free parameters of the models. The code has been recently made publicly available 1 .
Here we apply the method over a mock catalogue of dark matter haloes to forecast the constraints on the parameters in the two classes of models cited above, obtainable by current and future kinematic and lensing measurements of galaxy cluster mass profiles. This paper summarizes the results presented at the ALTECOSMOFUN'21 conference and is based upon the works published in 2020, 2021 and 2022.
In Section 2, we provide a brief overview of the MG models implemented, in Section 3 the main features of MG-MAMPOSST are presented while in Section 4 we show our results and preliminary outcomes from the application of MG-MAMPOSST over real data. Finally, a summary and a discussion about systematics are provided in Section 5

Theoretical Background
In scalar-tensor theories, the screening of the fifth force operates in different ways depending on the non-linear interactions of the field.
Models which implement the chameleon mechanism-typically conformally coupled models such as f (R) gravity, [14]-rely on the potential of the field to make the mass of the scalar very large at small scales, such that the fifth force does not propagate. In CS, the usual Newtonian potential is modified by an additional term which depends on the gradients of the scalar field φ (see e.g., [15,16]), where M(r) is the total mass enclosed within a sphere of radius r, M P = (8πG) −1/2 is the reduced Planck mass, G is the Newton's constant and Q is a dimensionless coupling constant 2 which can assume different values depending on the specific theory. The case of f (R) gravity, where the Einstein-Hilbert action is modified by adding a general non-linear function of the Ricci curvature scalar R, can be shown to be conformally equivalent to a subclass of chameleon models where Q = 1/ √ 6 (see e.g., [17]). In the above Equation (1) spherical symmetry has been assumed.
The explicit expression of the field profile depends on the matter density distribution ρ(r). Here we use the Navarro-Frenk-White (NFW) profile of [18] to model the matter density in MG frameworks, which has been shown to provide an overall good description of the total mass profile of cluster-size halos in equilibrium configuration from cosmological simulations and observational data in GR (e.g., [19]) and in MG (e.g., [20]). Nevertheless, the MG-MAMPOSST method can be easily extended to other mass models in non-standard frameworks. The NFW model is fully specified by the scale radius r s at which the logarithmic derivative of the profile is equal to −2 and the "virial" radius r 200 encloses an overdensity 200 times the critical density of the universe. Given the NFW model for ρ(r), we solve Equation (1) to obtain the field profile inside (where field gradients are negligilble) and outside the source (where the potential is subdominant) by assuming the same analytical approximation of e.g., [7], In the above equation ρ s (r s , r 200 ) is the characteristic density of the NFW model, φ ∞ is the background value of the field, C is an integration constant and x = r/r s ,. The solutions match at the screening radius S. Requiring the continuity of the function and of its first derivative at the matching point, C and S are fully specified in terms of (r 200 , r s , Q, φ ∞ ). In particular, the screening radius depends on r 3 200 . Thus, haloes with different masses exhibit a less/more efficient chameleon mechanism.
In CS, due to the conformal structure of the models, photon propagation is not affected by the fifth force. This can be seen by the definition of the lensing potential Φ lens which, at linear order, is given by the sum of the potential Φ and the relativistic potential Ψ. The contribution of the fifth force appears with opposite sign in the two terms and cancels out (see e.g., [16]). This means that lensing observations are sensitive only to the standard Newtonian potential.
In Vainsthein screening, which is a typical feature of DHOST theories, the suppression of the fifth force is achieved by working on higher-order derivatives of the field. In this case the mechanism to recover GR can be partly broken inside a massive object, giving rise to a fifth force which depends on the gradients of the matter distribution (see e.g., [21]). Assuming spherical symmetry and an NFW model for the matter density distribution, the Poisson equation associated with the Newtonian potential Φ (which governs the dynamics of galaxies and gas in clusters) is modified according to (e.g., [9]): where M NFW (r) is the NFW mass profile, M 200 is the mass of a sphere of radius r 200 enclosing an average density 200 times the critical density of the universe and c = r 200 /r s is called concentration. Finally, Y 1 is a dimensionless parameter describing the fifth force coupling. In this class of models, photon propagation is explicitly modified. In particular, the relativistic potential Ψ receives a contribution form the fifth force, giving rise to an effective lensing mass (e.g., [22]): where The coupling Y 2 appears only in the relativistic sector, i.e., it can be constrained only by lensing observations. Current constraints for Y 1 are of the order of 10 −2 as obtained from stellar probes (e.g., [23][24][25]), and of 10 −1 at the cosmological level with galaxy clusters (e.g., [26,27]). As for Y 2 , only cosmological bounds are available (O(1) from [9] and O(0.1) from [27], which worked on a generalization of the model presented here).

The MG-MAMPOSST Method
MG-MAMPOSST is a code aimed at constraining modified gravity models by analysing the kinematics of cluster member galaxies. Derived from the original MAMPOSST method 3 of [12], the program was first presented in a preliminary version in [11], and then extended and updated by [22]. The latest version implements both the CS and VS parametrization of the gravitational potentials described in Section 2, based on the NFW model for the matter density profile.
The input data-set of the MAMPOSST and MG-MAMPOSST procedure is the projected phase space (p.p.s.) of the member galaxies (R, v z ), where R is the projected radius from the cluster center and v z represents the line-of-sight (l.o.s.) velocity, computed in the rest frame of the cluster. Assuming dynamical relaxation (i.e. galaxies are collisonless tracers of the gravitational potential) and a Gaussian modelling of the three-dimensional velocity distribution, the codes computes orbits of cluster members by solving the (stationary) spherical Jeans' equation to obtain the velocity dispersion profile along the radial direction σ 2 r (see e.g., [28]): where ν(r) corresponds to the number density profile of tracers, β ≡ 1 − (σ 2 θ + σ 2 φ )/2σ 2 r is the velocity anisotropy profile and Φ is the total gravitational potential, which carries information about the nature of the gravitational interaction (e.g., Equations (1) and (3)).
Given a parametric expression for the above mentioned quantities, MG-MAMPOSST performs a Maximum Likelihood estimation of the model parameters using data of the p.p.s. In particular, the code computes the probability that a member galaxy found at the point (R i , v z,i ) in the p.p.s. belongs to the orbits distribution described by the model(s). In the most general case, the code can work with six free parameters, namely two mass profile parameters (e.g., r 200 , r s ), one parameter for the scaling of the number density profile r ν 4 , one parameter describing the velocity anisotropy profile β, and two parameters defining the MG model to be constrained (Q, φ ∞ for CS and Y 1 , Y 2 for VS).
The number density ν(r) can be, in general, excluded from the MG-MAMPOSST fit as it can be measured directly by analysing the projected distribution in the phase space. In the following, we assume an NFW profile to model the galaxy density profile. As for the velocity anisotropy, six possible parametrizations of β(r) are currently available in the code (see the technical manual [13]). As a case study, here we adopt the Tiret anisotropy model of [30], where β ∞ is the velocity anisotropy for r → ∞ and r β is the characteristic radius of β T (r) (anisotropy radius). In the current version of the code, r β is assumed to be equal to the scale radius of the mass profile r s . Moreover, hereafter we will work with the re-scaled anisotropy The parameter space can be explored either by computing the likelihood over a multidimensional grid of values, or by performing a Monte Carlo Markov Chain (MCMC) based on a simple Metropolis-Hastings algorithm with a fixed-step Gaussian random walk. The code takes few hours to produce a complete chain of ∼ 10 5 points.
Another feature introduced in MG-MAMPOSST is the possibility to simulate additional lensing information to be combined with the likelihood from internal kinematics. This is straightforward in the case of CS where lensing is not affected by the fifth force contribution; the lensing distribution is modeled as a Gaussian P lens (r s , r 200 ) for the NFW mass profile parameters r 200 , r s . The central values, the standard deviations and the correlation defining the distribution can be customized by the user. In the forecast analysis presented in Section 4, we will assume the lensing Gaussian is centered on the true values of the cluster mass profile's parameters.
The case of VS requires a full lensing simulation, as the lensing mass is explicitly modified in terms of the two fifth-force couplings Y 1 , Y 2 . In MG-MAMPOSST, we implemented a simple weak-lensing simulation which generates a mock reduced tangential shear profile assuming a fiducial NFW model in GR. The log-likelihood is given by where g t (R i ) is the simulated averaged reduced tangential shear profile at projected position R i , g t,vs (R i )|θ l ) is the theoretical profile computed in VS for the set of parameters θ l = (r s , r 200 , Y 1 , Y 2 ). The number of bins is fixed to N b = 10 (in agreement with current lensing surveys, e.g., [31]) in the projected radial range [0.12 R up , 2.9 R up ], where R up is the maximum value of the projected radius a galaxy can have to be considered in the MG-MAMPOSSt fit (see Section 4). The uncertainties are given by the quadratic sum of two contributions: where σ 2 e,i = σ 2 g /[π(α 2 up − α 2 low )n g ] is the noise due to the intrinsic ellipticity σ 2 g of the sources lying within an annulus between the angles α low and α up , and σ 2 lss is the uncertainty due to the projected large-scale structure. The average number of source galaxies per arcmin 2 is given by n g . The central values of the NFW parameters from which the shear profile is derived, as well as the lensing uncertainties and n g can manually set in the input file of MG-MAMPOSST.

Synthetic Halo Catalogue
The mock catalogue of dark matter haloes used to test MG-MAMPOSST has been produced with the CLUSTERGEN code (e.g., [32]), a generator of spherically symmetric, isolated distributions of particles in dynamical equilibrium, characterized by Gaussian 3D velocity distributions. We assume that all the systematics are under control, i.e., particles in each halo follow an NFW distribution and their velocity is assigned by assuming that σ 2 r is given by Equation (6). As for the other velocity dispersion components, We will comment on the effect of systematics in Section 5. It is important to point out here that a more rigorous way to compute orbits of a spherical system in dynamical equilibrium is to use six-dimensional distribution functions (e.g., [33,34]). A new version of CLUSTERGEN based on this methodology is currently under development.
We considered a total of 20 synthetic massive cluster-size haloes generated in GR 5 , as a reasonable number of relaxed galaxy clusters for which high-quality data could be available from upcoming surveys. All haloes are different realizations of the same NFW distribution with r 200 = 2.0 Mpc and r s = 0.3 Mpc. For the case of CS, since the true values of the NFW parameters affect the results, we will further discuss how the constraints change when varying the mass of the synthetic cluster.
For each halo we generated two p.p.s., obtained by considering 600 and 100 particles within the radial range [0.05 Mpc, r 200 ], as an optimistic expectation of the number of member galaxies' spectroscopic redshifts available from upcoming surveys, although not unrealistic. The bounds of the radial range are set in spite of real observations, to exclude the cluster core where the Brightest Central Galaxy (BCG) dominates the internal dynamics. As for the upper bound, we adopt the conservative limit of R up = r 200 to ensure the validity of the Jeans' equation. Finally, as mentioned in Section 3, we consider a Tiret model for the velocity anisotropy profile with β ∞ = 0.5 (i.e. A ∞ = 1.41).

Vainsthein Screening
We apply the MG-MAMPOSST method to the synthetic p.p.s. with additional simulated lensing information, assuming the VS parametrization for the gravitational potentials, Equation (3) and Equation (4), to constrain the set of parameters r 200 , r s , A ∞ , Y 1 , Y 2 . Assuming uniform uninformative priors on each parameter, we perform an MCMC sampling of the joint likelihood: is the total MG-MAMPOSST likelihood obtained by combining the information of N h phase spaces, with N h = 1 . . . 20.
As for the lensing simulation, we have set n g = 30 arcmin −2 , as expected for the Wide Survey of the Euclid mission [35], σ lss = 0.005 and σ g = 0.3.
The results of our forecast are summarized in Figure 1 for one p.p.s. with 600 tracers and in columns two to four of Table 1. While for Y 1 we cannot obtain competitive bounds with respect to stellar probes, even when increasing the number of haloes in the fit, the relativistic coupling Y 2 can be constrained at the level of O(0.1) already with a few clusters, in agreement with the results of [27] and with a ∼ 2 − 3 times improvement with respect to the analysis of [9]. In particular, for a single halo we obtain Y 2 = 0.08 +0. 32 −0.28 (600 tracers) and Y 2 = 0.10 +0.44  The GR expectation (i.e., Y 1 = Y 2 = 0) is indicated by the black vertical dashed lines. From [22]. Table 1. Constraints at 95% C.L obtained for the parameters of the two modified gravity models presented in this work applying the MG-MAMPOSST fit with additional lensing information over the mock catalogue of synthetic halos with r 200 = 2.0 Mpc and r s = 0.3 Mpc. From column two to column six: Vainshtein screening. Column seven and eight: scalaron field in general chameleon f (R) gravity, related to chameleon field φ ∞ through Equation (13). For both models we show the results when using ∼600 and ∼100 cluster members in the MG-MAMPOSST fit. As a preliminary application of MG-MAMPOSST over real data, we further use the highprecision kinematic and lensing information for the massive relaxed cluster MACS J1206.2-0847 (MACS 1206 hereafter) at redshift z = 0.44, analysed within the Cluster Lensing and Supernova Survey with Hubble (CLASH, [36]) and CLASH-VLT [37] collaborations, to constrain Y 1 and Y 2 in VS (see [38] for details). In this case, we have combined the Jeans' analysis of 375 cluster member galaxies in the p.p.s. with the strong plus weak lensing data of [31]. The results of

Vainshtein Screening f (R) Gravity
are in very good agreement with the prediction of our forecast analysis. In the above equation, the systematic uncertainties encapsulate the effect of changing the parametrization of the velocity anisotropy profile and the value of the scale radius of the number density profile r ν which, as mentioned above, is given by an external fit of the p.p.s.

Chameleon Screening
In this case, the MG parameters defining the model are φ ∞ and Q. As for VS, we perform an MCMC sampling of the joint kinematic and lensing likelihood. In the Gaussian distribution P lens (r 200 , r s ), we assume average realistic uncertainties for current lensing surveys σ r 200 /r 200 = 0.1 and σ r s /r s = 0.3, with a correlation ρ = 0.5. As in e.g. [7,8], we work with the rescaled variables which spawn the range [0, 1]. In Figure 2 we have plotted the resulting 2σ (darker areas) and 3σ (lighter areas) allowed regions in the space (φ 2 , Q 2 ) for three relevant cases. These results are very similar to that of [7,8] who performed joint X-rays and weak lensing analyses of galaxy clusters, but in the case of MG-MAMPOSST plus lensing, the excluded region is shifted towards larger values of Q 2 . This difference is related to the fact that hot X-ray-emitting gases in clusters and member galaxies are subject to distinct physical processes, despite the fact that they perceive the same gravitational potential. Note that the constraints show only a mild dependence on the number of tracers in the p.p.s., as the strength of the results relies on the additional (lensing) information on r s and r 200 , needed to break the statistical degeneracy between the mass profile parameters and the MG parameters in internal kinematic analysis. When increasing the number of haloes in the fit (left panel in Fig. 2, the bounds on the allowed region are noticeably tightened. As a last step, we focus on the popular f (R) gravity subclass of chameleon models, where Q = 1/ √ 6. In these framework, the fifth force is carried out by the derivative of the function with respect to the Ricci scalar, f R = ∂ f /∂R, which acts as a scalar field generally known as scalaron. The present-day background value of the scalaron f R0 is connected to the chameleon field φ ∞ through (e.g., [8,17]): where c l is the speed of light. We perform the same kinematic and lensing analysis as for the general chameleon case described above, to constrain the value of the background scalaron field. It is worth pointing out that the constraints on f (R) gravity derived with this approach are independent of the choice of the function f (see [7,22]). The marginalized distributions of φ ∞ obtained in different cases are shown in Figure 13 of [22] for the combined analysis of 10 haloes. We further list the constraints for f R0 at 95% C.L. in the last two columns of Table 1 when varying the number of haloes considered in the fit. As in CS framework, the screening radius depends on the mass and the density of the sources (i.e., on r s and r 200 ), we have also generated p.p.s. from a less massive halo with a smaller concentration, characterized by r 200 = 1.41 Mpc and r s = 0.33 Mpc. The marginalized 10-halo distribution is shown by the blue curve in Figure 13 of citepPizzuti2021 What we found is that, as expected, smaller clusters exhibit a less efficient screening mechanism, allowing us to better constrain possible departures from GR.

Discussion and Conclusions
We have introduced the main features of the license-free code MG-MAMPOSST, developed to constrain modified gravity models with the internal kinematic analysis of galaxy clusters, further equipped with simple lensing simulations for the two classes of scalar tensor theories currently implemented, chameleon screening and Vainsthein screening. The tests performed over mock haloes show that the combination of internal kinematic and lensing analyses provides promising results, even for a limited number of clusters. In particular, the relativistic coupling Y 2 in Vainsthein screening can be constrained with uncertainties of a few percent with a dozen clusters, assuming reliable high-quality kinematics and lensing information. The exercise has been pushed forward by using the currently available data of the massive galaxy cluster MACS 1206; in this case, we obtain Y 2 = −0.12 +0.66 −0.67 (see [38]), in agreement with GR prediction and comparable to other constraints at cluster scales for these models [9,27].
In a similar fashion, the allowed region in the parameter space of chameleon gravity (φ 2 , Q 2 ) can be strongly tightened with joint lensing and internal kinematic analyses of a few galaxy clusters, in a way which is complementary to cluster mass determinations using Xray analyses of the hot intra-cluster gas (e.g., [7,8]). This indicates that the combination of kinematics of member galaxies, X-rays and lensing probes could provide very strong constraints on the behavior of gravity at cluster scales.
So far, we have investigated the constraining power of the MG-MAMPOSST method, assuming that all the systematics are under control. However, deviations from the dynamical relaxation assumption, crucial for the Jeans' analysis in clusters, as well as departures from spherical symmetry, can in principle affect the bounds on the model parameters, giving rise to spurious detection of modified gravity. Calibration and estimation of systematic effects is fundamental in view of the large amount of imaging and spectroscopic data that will be available in the near future. In [39], we showed that the analysis of cluster-size dark matter haloes in a ΛCDM universe indicates departures from GR in ∼70% of cases. Nevertheless, we also found that this percentage can be drastically decreased by looking at the distribution of member galaxies in the p.p.s and the line-of-sight velocity distribution P(v z ). In particular, deviations from Gaussianity in P(v z ), quantified by the so-called Anderson-Darling coefficient A 2 [40], are strongly related to the probability of finding suitable clusters for the application of our method. This is the case of MACS 1206, characterized by a small value of A 2 ∼ 0.6, in agreement with a nearly Gaussian velocity distribution.
Given the current capabilities and limitations of our method, several improvements can be made, which will be available within the next versions of the MG-MAMPOSST code. In particular, gravitationally bounded structures in some modified gravity theories may not be correctly described by an NFW profile (see, e.g., [41] and references therein). As such, we are planning to extend the parametrizations of the modified gravitational potentials to other mass models which may provide a better description of haloes in MG frameworks. Moreover, the efficiency of the screening mechanism depends on the shape of the matter distribution (e.g., [42]); thus, one should go beyond the assumption of spherical symmetry in computing the orbits with the Jeans' equation, in order to fully characterize the dynamics in these MG models. Finally, our forecast analysis has been performed on synthetic phase spaces generated in GR to explore the constraining power of MG-MAMPOSST. More realistic constraints can be obtained by computing orbits of member galaxies directly in a modified gravity scenario and comparing the results to real observational data. This can be conducted with the new version of the CLUSTERGEN code we are developing, which will be made publicly available in the future together with MG-MAMPOSST.
Funding: LP is partially supported by a 2019 "Research and Education" grant from Fondazione CRT. The OAVdA is managed by the Fondazione Cleḿent Fillietroz-ONLUS, which is supported by the Regional Government of the Aosta Valley, the Town Municipality of Nus and the Unitedes Communes valdotaines Mont-Eḿilius.

Informed Consent Statement: Not applicable
Data Availability Statement: The MG-MAMPOSST code is publicly available. A test p.p.s. data-set from a synthetic halo is shipped with the code. The data of MACS 1206 are partially provided by [31,43] by permission. Data will be shared on reasonable request to the corresponding author with the permission of CLASH-VLT collaboration: [37]. Notes 1 https://github.com/Pizzuti92/MG-MAMPOSSt. Accessed 19 February 2022, see [13] for the basic usage of the code. 2 In the literature the coupling constant is often indicated by β. However, as β(r) also denotes the velocity anisotropy profile in kinematic analyses of galaxy clusters, we adopt Q for the coupling to avoid confusion. 3 The public version of MAMPOSST can be found at https://gitlab.com/gmamon/MAMPOSSt, accessed 19 February 2022. 4 Note that r ν is in general different from r s as the distribution of galaxies in clusters may not follow the distribution of the total matter (e.g., [29]). 5 In principle CLUSTERGEN can be used to produce mock clusters adopting different modified gravity setups and matter density distributions. In the exercise presented here we focus only on an NFW profile in GR as our fiducial model.