Future constraints on the gravitational slip with the mass profiles of galaxy clusters

TL;DR

This study demonstrates that galaxy cluster mass profiles from lensing and galaxy dynamics can constrain the gravitational slip parameter to percent-level accuracy with current and future surveys, aiding tests of gravity theories.

Contribution

It introduces a method combining lensing and dynamical data to forecast constraints on the gravitational slip parameter using galaxy clusters.

Findings

Galaxy clusters can constrain the gravitational slip to percent accuracy.

Current and future surveys can achieve these constraints with tens of clusters.

Systematic effects have a mild impact on the constraints.

Abstract

The gravitational slip parameter is an important discriminator between large classes of gravity theories at cosmological and astrophysical scales. In this work we use a combination of simulated information of galaxy cluster mass profiles, inferred by Strong+Weak lensing analyses and by the study of the dynamics of the cluster member galaxies, to reconstruct the gravitational slip parameter $\eta$ and predict the accuracy with which it can be constrained with current and future galaxy cluster surveys. Performing a full-likelihood statistical analysis, we show that galaxy cluster observations can constrain $\eta$ down to the percent level already with a few tens of clusters. We discuss the significance of possible systematics, and show that the cluster masses and numbers of galaxy members used to reconstruct the dynamics mass profile have a mild effect on the predicted constraints.

Tables

Table 1. Table 1 : A summary with the characteristics of the synthetic clusters sample used in our analysis ordered according to increasing mass. The second column shows the parameters r 200 subscript 𝑟 200 r_{200} for each cluster, sampled with a fixed step size of 0.15 ​ Mpc 0.15 Mpc 0.15\,\mbox{\rm Mpc} , while the third column the corresponding r s subscript 𝑟 𝑠 r_{s} parameter derived according to the empirical concentration-mass relation ( 4.3 ). The fourth column shows the scale radius r ν subscript 𝑟 𝜈 r_{\nu} , obtained by fitting the projected number density profile from the phase space with a projected NFW model, which is additionally required in the MAMPOSSt procedure. The fifth column denotes the number of tracers within a sphere of radius r = r 200 𝑟 subscript 𝑟 200 r=r_{200} , computed according to eq. ( 4.2 ).

Cluster ID	${𝐫}_{\mathrm{𝟐𝟎𝟎}}\left[\text{Mpc}\right]$	${𝐫}_{𝐬}\left[\text{Mpc}\right]$	${𝐫}_{\nu }\left[\text{Mpc}\right]$	# of tracers ${𝐍}_{𝐠}\left({𝐫}_{\mathrm{𝟐𝟎𝟎}}\right)$
1	0.750	0.185	0.169	300
2	0.900	0.094	0.097	324
3	1.050	0.177	0.178	343
4	1.200	0.133	0.144	384
5	1.350	0.231	0.219	437
6	1.500	0.258	0.256	500
7	1.650	0.227	0.252	570
8	1.800	0.360	0.387	648
9	1.950	0.511	0.466	732
10	2.100	0.515	0.473	823
11	2.250	0.539	0.503	920
12	2.400	0.614	0.641	1024
13	2.550	0.848	0.807	1133
14	2.700	0.993	1.005	1249
15	2.850	1.112	1.117	1371

Table 2. Table 2 : Constraints from the combined analysis of multiple number of clusters for scale-independent (third column) and scale-dependent (fourth column) gravitational slip η 𝜂 \eta . In the latter case, the error is evaluated at the reference scale r = 1.5 ​ Mpc 𝑟 1.5 Mpc r=1.5\,\text{Mpc} . In each block each of the two rows refer to a different sub-case: the first one shows the constraints from the 15 clusters-likelihood, while the second refers to the extrapolated constraints for N c ​ l ​ u ​ s ​ t ​ e ​ r = 75 subscript 𝑁 𝑐 𝑙 𝑢 𝑠 𝑡 𝑒 𝑟 75 N_{cluster}=75 . When varying the lensing parameters (correlation ρ 𝜌 \rho and relative errors respectively), the number of tracers in the dynamics fit is fixed to N g = 1000 subscript 𝑁 𝑔 1000 N_{g}=1000 for all clusters.

case	$N_{cluster}$	$\eta =\text{constant}$	$\eta =\eta \left(r=1.5\text{Mpc}\right)$
reference	15	$1.00\pm 0.05\left(1\sigma \right)\pm 0.09\left(2\sigma \right)$	${1.00}_{-0.09}^{+0.10}\left(1\sigma \right)\pm 0.21\left(2\sigma \right)$
	75	$1.00\pm 0.02\left(1\sigma \right)\pm 0.04\left(2\sigma \right)$	${1.00}_{-0.03}^{+0.04}\left(1\sigma \right)\pm 0.08\left(2\sigma \right)$
Variable $N_g$	15	$1.00\pm 0.05\left(1\sigma \right)\pm 0.09\left(2\sigma \right)$	$1.00\pm 0.12\left(1\sigma \right){}_{-0.24}{}^{+0.25}\left(2\sigma \right)$
	75	$1.00\pm 0.02\left(1\sigma \right)\pm 0.04\left(2\sigma \right)$	$1.00\pm 0.04\left(1\sigma \right){}_{-0.09}{}^{+0.10}\left(2\sigma \right)$
$\rho =0.15$	15	$1.00\pm 0.06\left(1\sigma \right),\pm 0.11\left(2\sigma \right)$	$1.00\pm 0.13\left(1\sigma \right),\pm 0.27\left(2\sigma \right)$
	75	$1.00\pm 0.03\left(1\sigma \right),\pm 0.05\left(2\sigma \right)$	$1.00\pm 0.05\left(1\sigma \right)\pm 0.09\left(2\sigma \right)$
$\rho =0$	15	${1.00}_{-0.05}^{+0.06}\left(1\sigma \right),\pm 0.10\left(2\sigma \right)$	$1.00\pm 0.12\left(1\sigma \right),\pm 0.26\left(2\sigma \right)$
	75	${1.00}_{-0.02}^{+0.03}\left(1\sigma \right),\pm 0.04\left(2\sigma \right)$	$1.00\pm 0.05\left(1\sigma \right)\pm 0.08\left(2\sigma \right)$
$\rho =0.99$	15	$1.00\pm 0.03\left(1\sigma \right),\pm 0.11\left(2\sigma \right)$	$1.00\pm 0.08\left(1\sigma \right),\pm 0.19\left(2\sigma \right)$
	75	$1.00\pm 0.01\left(1\sigma \right),\pm 0.03\left(2\sigma \right)$	$1.00\pm 0.03\left(1\sigma \right)\pm 0.06\left(2\sigma \right)$
$\sigma \left(r_{200}\right)=0.15\times r_{200}$	15	${1.00}_{-0.10}^{+0.11}\left(1\sigma \right){}_{-0.17}{}^{+0.18}\left(2\sigma \right)$	$1.00\pm 0.15\left(1\sigma \right)\pm 0.30\left(2\sigma \right)$
	75	$1.00\pm 0.05\left(1\sigma \right)\pm 0.08\left(2\sigma \right)$	$1.00\pm 0.05\left(1\sigma \right)\pm 0.11\left(2\sigma \right)$
$\sigma \left(r_s\right)=0.45\times r_s$	15	$1.00\pm 0.05\left(1\sigma \right)\pm 0.09\left(2\sigma \right)$	$1.00\pm 0.15\left(1\sigma \right)\pm 0.31\left(2\sigma \right)$
	75	$1.00\pm 0.02\left(1\sigma \right)\pm 0.04\left(2\sigma \right)$	$1.00\pm 0.05\left(1\sigma \right)\pm 0.12\left(2\sigma \right)$

Equations

$$\frac{\partial (\nu {\sigma }_r^2)}{\partial t}+2\beta (r)\frac{\nu {\sigma }_r^2}{r}=-\nu (r)\frac{\partial \Psi }{\partial r},$$

$${\nabla }^2\Psi =4\pi G{\rho }_m,$$

$$M_{\text{tot}}(R)=4\pi {\int }^R{r}^2{\rho }_m(r)\text{mbox}dr.$$

$${\rho }_{\text{NFW}}(r)=\frac{{\rho }_{\text{s}}}{(r/r_s)(1+r/r_s{)}^2},$$

$${\rho }_s={\rho }_c(z)\frac{{\Delta }_{vir}}{3}\frac{c^3}{\log (1+c)-c/(1+c)}.$$

$$M(r)=M_{200}\frac{\ln (1+r/r_s)-(r/r_s)(1+r/r_s{)}^{-1}}{\ln (1+c)-c/(1+c)},$$

$$\Psi (r)=G{\int }_{r_0}^r\frac{ds}{s^2}{M}_{\text{dyn}}(s).$$

$${\nabla }^2(\Phi +\Psi )=8\pi G{\rho }_{\text{lens}},$$

$$M_{\text{lens}}=\frac{r^2}{2G}\frac{d}{dr}(\Phi +\Psi ).$$

$$\Phi (r)=G{\int }^r\frac{ds}{s^2}\left[2M_{\text{lens}}(s)-M_{\text{dyn}}(s)\right].$$

$$\eta (r)=\frac{{\int }^r\frac{ds}{s^2}\left[2M_{\text{lens}}(s)-M_{\text{dyn}}(s)\right]}{{\int }^r\frac{ds}{s^2}{M}_{\text{dyn}}(s)}.$$

$$\eta =h_2\cdot \frac{1+h_4{k}^2}{1+h_5{k}^2},$$

$$\eta (r,t)=f_1(t)\left(1+f_2(r,t)\right),$$

$$\eta (r;r_{200}^D,r_s^D)\ne \eta (r;r_{200}^L,r_s^L),$$

$$P_L(r_s,r_{200})=\frac{1}{2\pi {\sigma }_{r_s}{\sigma }_{r_{200}}\sqrt{1-{\rho }^2}}\exp \{-\frac{1}{2(1-{\rho }^2)}[\frac{(r_s-{\bar{r}}_s{)}^2}{{\sigma }_{r_s}^2}+$$

$$+\frac{(r_{200}-{\bar{r}}_{200}{)}^2}{{\sigma }_{r_{200}}^2}-\frac{2\rho (r_s-{\bar{r}}_s)(r_{200}-{\bar{r}}_{200})}{{\sigma }_{r_s}{\sigma }_{r_{200}}}]\},$$

$$N_{g}^{(i)}(r_{200})=\left\{\begin{array}[]{lll}300&&\;\text{$i=1$}\,,\\ \left(M^{(i)}_{200}/M^{(1)}_{200}\right)\cdot\frac{N_{c}}{i}&&\;\text{$i>1$}\,.\end{array}\right.$$

$$c(M_{200},z)=A\times {\left(\frac{1.37}{1+z}\right)}^B\cdot {\left(\frac{M_{200}}{8\times 10^{14}{M}_{\odot }/h}\right)}^K,$$

$$A=3.66\pm 0.16,B=-0.14\pm 0.52,K=-0.32\pm 0.18.$$

$$\beta (r{)}_T={\beta }_{\infty }\cdot \frac{1}{1+r_c/r},$$

$${\sigma }_r^2(r)=\frac{1}{\nu (r)}{\int }_0^{\infty }\exp \left[2{\int }_r^s\beta (t)\frac{dt}{t}\right]\nu (s)\frac{d(\Psi (s))}{ds}ds,$$

$${\sigma }_{\theta }^2(r)=[1-\beta (r)]{\sigma }_r^2(r).$$

$$q(R,v_z)=\frac{2\pi Rg(R,v_z)}{N_p(R_{max})-N_p(R_{min})}.$$

$$g(R,v_z)=\sqrt{\frac{2}{\pi }}{\int }_R^{\infty }\frac{r\nu (r)}{\sqrt{r^2-R^2}}\frac{dr}{{\sigma }_r(r)\sqrt{1-\beta (r)R^2/r^2}}\exp \left[-\frac{v_z^2}{2(1-\beta (r)R^2/r^2){\sigma }_r^2(r)}\right].$$

$$-\log \mathcal{L}=-\sum _{i=1}^n\log q(R_i,v_{z,i}|\vec{\theta }),$$

$$\eta =2\times \frac{\ln (1+c^{\text{dyn}})-c^{\text{dyn}}/(1+c^{\text{dyn}})}{\ln (1+c^{\text{lens}})-c^{\text{lens}}/(1+c^{\text{lens}})}{\left(\frac{r_{200}^{\text{lens}}}{r_{200}^{\text{dyn}}}\right)}^3-1,$$

$$\eta =1.00\pm 0.05\text{at}1\sigma ,$$

$$\eta (N_g=100)=1.00\pm 0.07(1\sigma ){}_{-0.11}^{+0.13}(2\sigma )\text{15 clusters},$$

$$\eta (N_g=100)=1.00\pm 0.03(1\sigma ){}_{-0.05}^{+0.06}(2\sigma )\text{75 clusters}.$$

Full text

1 Introduction

The dark energy problem has challenged our understanding of gravity at large scales in the Universe. To explain it, a variety of different theoretical scenarios have been proposed based on the introduction of new degrees of freedom beyond the standard paradigm of General Relativity (GR) (for a review see e.g [1, 2]). The recent impressive measurement of the speed of gravitational waves provided the most stringent test on viable modified gravity theories so far, with direct implications for cosmology and the large scale structure of the Universe. However, understanding the nature of the late–time acceleration of the Universe calls for further, complementary observational information, that could shed light on the viability of the surviving theory space.

From a phenomenological point of view, modifications of GR leave two key imprints at the level of large scale structures in the Universe. In particular, they affect the effective strength of gravity, impacting the clustering of matter at large scales, while they may also impact the way gravity interacts with light, leaving a characteristic imprint on weak lensing observables. Both effects, however, could be equally due to new gravitational degree(s) of freedom that do not modify gravity, such as minimally-coupled scalar-tensor models (quintessence, k–essence), so that they cannot be used to detect genuine departures from GR. A crucial parameter in this regard is the gravitational slip $\eta$ defined as the ratio between the effective gravitational coupling of light to that of matter (to be defined later on). In GR and in models where the new degree of freedom is minimally-coupled to gravity one has $\eta=1$, while an $\eta\neq 1$ would be the smoking gun for a departure from GR (see e.g [3]) at late times, at least as long as dark matter can be described as a perfect fluid 111We notice here that, free streaming neutrinos can also be a source of gravitational slip. However, at late times in the cosmological evolution their contribution can be neglected.. Crucially, $\eta$ can be reconstructed in a model-independent manner by combining galaxy kinematics and weak lensing observations [4, 5], while it is directly linked to the properties of gravitational wave propagation at any scale [6, 7, 8]. Therefore, its reconstruction from observations is key in testing departures from GR and/or distinguishing amongst large classes of gravity models at both cosmological and astrophysical scales.

Galaxy clusters are ideal laboratories in this regard, as they can provide us with simultaneous information about the local gravitational potentials sourcing the dynamics of their member galaxies as well as the lensing of light, which can be combined to reconstruct $\eta$. In this sense, tests of gravity based on $\eta$ can be understood as a consistency check between the cluster’s mass profile as inferred from kinematical and lensing observations respectively. Current and future galaxy cluster surveys observing at different, complementary bands, such as Euclid, Dark Energy Survey (DES), XMM-Newton, or eROSITA 222http://www.darkenergysurvey.org/, http://www.mpe.mpg.de/heg/www/Projects/erosita/index.php, http://sci.esa.int/xmm-newton/43133-slew-survey-and-catalogue/. For a review on future surveys see [9]. aim to provide observations of hundreds of clusters with accurate mass determination. This work adds to previous efforts to constrain gravity with galaxy clusters in various theoretical setups, e.g in [10, 11, 12, 13], and it is complementary to conceptually similar reconstructions of $\eta$ using redshift space distortions, galaxy-galaxy lensing and the background expansion [14].

The goal of the present work is two-fold: On the one hand, we aim at predicting the potential statistical accuracy of future constraints on $\eta$ from galaxy cluster kinematics and lensing, while at the same time understanding their dependence on key factors such as the number of member galaxies and the mass of the clusters in the sample. Our analysis is based on the dynamical mass profile reconstructions, performed with the MAMPOSSt method of Ref. [15], of a set of simulated spherically symmetric isolated clusters. This approach will give us enough flexibility to investigate how different cluster masses and densities of tracers affect our constraints. For the complementary lensing information we will rely on the lensing mass profile inferred by Refs. [16, 17, 18] on the massive, relaxed galaxy cluster MACS 1206, which belongs to a sample of 20 X-ray selected objects analysed within the CLASH/ CLASH-VLT collaborations (refs. [19, 20]).

It should be noted that theories modifying gravity typically introduce a characteristic time and scale dependence on the gravitational slip parameter $\eta=\eta(t,k)$, which is often captured through appropriate parametrisations. Unlike previous constraints of gravity [10, 11, 12, 13] based on particular gravity setups, in our analysis we will make no assumption on the particular functional dependence of $\eta$ on time and scale, the only assumption being that of a Navarro-Frenk-White parametrisation for the cluster mass profile. In this regard, our constraints will be sufficiently general – any departures from the expectation $\eta=1$ is incorporated in a difference between the mass profile parameters inferred by dynamics and lensing probes, which can then be translated into constraints on the parameter space of particular theories of gravity.

We show that already with $\sim 15$ clusters we can get at $1\sigma$ an accuracy of $\leavevmode\nobreak\ \sim 5\%$ for a scale-independent and $\leavevmode\nobreak\ \sim 20\%$ for a scale-dependent $\eta$. The accuracy is significantly improved down to the few percent level for both cases when extrapolated to $\sim 75$ clusters. We investigate the effect of varying the number of tracers (i.e. member galaxies) in the dynamics analysis, the dependence of constraints on the cluster mass, as well as the uncertainties in the simulated lensing probability distribution.

We structure the paper as follows: Section 2 lays down the main assumptions and definitions of our analysis, while in Section 3 we introduce the essential theoretical framework. A detailed exposition of all our simulations is presented in Section 4, and finally, the results of the statistical analysis are explained in Section 5. The main results are summarised in the last section of the paper, while a summary of the predicted constraints can be found in Table 2.

2 Essential assumptions and definitions

An essential step towards forecasting the potential of future galaxy-cluster surveys to constrain the gravitational slip parameter is having a handle upon a measurement of clusters dynamical and lensing masses, as we explain in more detail in Section 3. In this regard, we will combine information from a set of reliable simulated data based on current and future observations. The analysis will inevitably rely on certain fundamental simplifications and assumptions about the structure and dynamics of the clusters and the underlying cosmology, which is important to lay out before we proceed. They can be summarized as follows:

i. The clusters will be modelled as spherical inhomogeneities in a complete equilibrium, that is, in thermal and hydrostatic equilibrium. In other words, the clusters will be assumed to be relaxed, with their mass distribution and gravitational potentials being time-independent. Although from an observational viewpoint it might be challenging to identify truly relaxed clusters, this simplification will allow us to avoid the introduction of systematics and/or degeneracies that we would otherwise encounter when dealing with a non-relaxed cluster. We comment on the issue of departures from relaxation in Section 3. Additional assumptions related to the anisotropy of the cluster’s velocity field are also explained in Section 4.1.

ii. The density distribution of the cluster will be dominated by dark matter, which we will model by a Navarro-Frenk-White (NFW) profile. The member galaxies will be assumed to be collissionless and non-interacting point particles moving along geodesics and populated according to the same profile. The validity of this assumption is ensured within the linear regime, whereas it is expected to break down as soon as non-linear galactic interactions become important. These are however rare, given the high speed of galaxy-galaxy encounters in clusters.

iii. The background cosmology will be that of the $\Lambda$CDM model, while large-scale structures will be sourced by linear scalar fluctuations according to $ds^{2}=-(1+\Psi(t,{\bf x}))dt^{2}+a(t)^{2}(1+\Phi(t,{\bf x}))d{\bf x}^{2}$. Under this convention, we define the gravitational slip parameter as $\eta\equiv\Phi/\Psi$.

iv. The functional form of the gravitational slip parameter $\eta=\eta(t,k)$ as a function of time and scale depends on the particular gravity action under consideration. Here, we shall make no a priori assumption about the actual gravity model, proceeding in a sufficiently general manner considering two distinct cases: A case where $\eta$ is assumed to be scale-independent throughout our clusters sample, and another one where it depends on the cluster’s spatial scale through the characteristic scale $r_{s}$ of the NFW profile function. These two cases encompass the predictions of scalar–tensor theories at different (subhorizon) scales, as explained in Section 3.

3 Gravitational slip as a consistency condition on the cluster’s dynamical and lensing mass

Testing the gravitational slip at the galaxy cluster level can be understood as a consistency check between the lensing and dynamical mass of the cluster, $M_{dyn}$ and $M_{lens}$ respectively. In GR, $M_{dyn}=M_{lens}$, which is no longer true as soon as gravity is modified 333Our definition of modified gravity is that of [8], defined as any theory modifying the dynamics of the genuine degrees of freedom of GR, that is the two polarisations of the spin-two field..To understand why this is so, let us assume a cluster to be in a state of dynamical equilibrium with the total gravitational potential $\Psi$. If interactions and dissipative effects are negligible, the cluster member galaxies acts like collisionless tracers of the gravitational potential and their dynamics is governed by the Jeans equation:

[TABLE]

with $\nu(r)$ the number density of tracers, $\sigma^{2}_{r}$ the velocity dispersion along the radial direction and $\beta\equiv 1-(\sigma_{\theta}^{2}+\sigma^{2}_{\phi})/2\sigma^{2}_{r}$ the velocity anisotropy profile. The gravitational potential $\Psi$ is connected to the matter density distribution through the Poisson equation:

[TABLE]

where $\rho_{m}$ stands for the total mass density of the cluster at a given radius $r$. The total density is assumed to be dominated by the dark matter component, i.e $\rho_{m}=\rho_{\text{DM}}+\rho_{\text{gas}}+\rho_{\text{gal}}\simeq\rho_{\text{DM}}$, as baryons provide a negligible contribution. The total mass of the (spherical) cluster with total radius $R$ is given by the mass conservation equation:

[TABLE]

In order to model the dark matter density we will follow a parametric approach through the use of a NFW profile [21], which describes the distribution of collisionless Cold Dark Matter (CDM) in spherical cluster-size halos in an equilibrium configuration, defined through

[TABLE]

with $\rho_{s}$ and $r_{s}$ two free parameters to be fit from observations. The NFW profile is in very good agreement with N-body simulations. In particular, $r_{s}$ is a characteristic scale radius, corresponding to the distance from the center of the cluster where $d(\ln\rho)/d\ln r=-2$, while $\rho_{s}$ is defined as

[TABLE]

Here, $\rho_{c}(z)=3H^{2}(z)/8\pi G$ is the critical density of the Universe at redshift $z$, and $c\equiv r_{vir}/r_{s}$ is the so–called concentration parameter, while $r_{vir}$ denotes the radius enclosing a mean overdensity $\Delta_{vir}$ times the critical density of the Universe. In this paper we will use $\Delta_{vir}=200$ (and so, $r_{vir}=r_{200}$), which is close to the density at virialization predicted by the spherical collapse model in a EdS universe $\Delta_{sc}=178$.

By inserting equation (3.5) into (3.3) we obtain the halo mass as a function of radius as

[TABLE]

with $M(r_{200})=M_{200}$ the total mass enclosed in a sphere of radius $r_{200}$. Despite its relative simplicity, the NFW profile has been widely adopted in the literature to describe the mass profile of gravitationally bound structures, as it has been shown to provide a generally good fit to observational data (e.g. refs. [22, 17]). The choice of a parametric description of the density perturbations has the noticeable advantage to simplify the extension of the analysis to specific modified gravity theories, where departures from GR can be easily incorporated in additional parameters in the model, as we will discuss in more detail later on. Now, assuming spherical symmetry the Poisson equation can be integrated to yield

[TABLE]

The above equation provides a definition for a cluster dynamical mass, that is, the mass one would infer through galaxy kinematics measurements, related to the Newtonian potential $\Psi$. We can derive a similar definition for the lensing mass, considering that light geodesics respond to the weak lensing potential $\Phi+\Psi$, described by the Poisson-type equation through

[TABLE]

which under the assumption of spherical symmetry becomes

[TABLE]

Equations (3.8) and (3.7) allow us to solve for the potential $\Phi$,

[TABLE]

Combining the expressions for the potentials $\Phi$ and $\Psi$ we can compute the ratio $\Phi/\Psi$, to arrive to a relation for the gravitational slip parameter $\eta$ in terms of the lensing and dynamical mass of the cluster as,

[TABLE]

Equation (3.11) will be our starting point in predicting the constraints on $\eta$ at the cluster level. Notice that, in GR it is $M_{\text{dyn}}=M_{\text{lens}}$ and $\eta=1$.

In the usual description of $\eta$ in general classes of gravity models, one usually works in momentum (Fourier) space, where the scale dependence is controlled through the Fourier wavenumber $k$. For example, for the popular class of the general Horndeski scalar-tensor theories in the sub-horizon quasi-static regime $\eta$ acquires the form

[TABLE]

with the functions $h_{i}=h_{i}(t)$ defined in [23]. For the particular case of Brans–Dicke/$f(R)$ gravity they depend on the scalaron Compton wavelength. However, since our simulations, and hence our analysis, will be in real space it is essential to discuss how our parametrisation relates to the Fourier space analyses. To get the real space $\eta$ one would have to express the linearised equations for $\Phi$ and $\Phi+\Psi$ in real space, assuming a parametrisation for the density profile. In general, we expect the resulting expression for for generic theories of gravity to acquire the following schematic form

[TABLE]

with GR corresponding to the limit $f_{1}\to 1,f_{2}\to 0$. A scale–independent $\eta$ is sourced only by $f_{1}$, which is in principle time-dependent – in our analysis, we will confine ourselves to a fixed redshift, which enforces $\eta$ to be a constant for this case. Without the assumption of a model-dependent parametrisation our approach will be rather phenomenological, yet sufficiently general to encompass any gravity theory. In particular, in view of (3.11) and the assumption of a NFW profile for the matter distribution, it is easy to see that $\eta$ will be sourced through a difference of the inferred parameters $r_{200},r_{s}$ from lensing and dynamics respectively, i.e

[TABLE]

unless $(r_{200}^{D},r_{s}^{D})=(r_{200}^{L},r_{s}^{L})$. From an observational viewpoint, the existence of some new degree of freedom sourcing $\eta\neq 1$, will manifest itself as a tension between the inferred NFW parameters from lensing and dynamics respectively, when trying to reconstruct the mass distribution of the cluster. In this regard, the strategy in our analysis will be as follows: the synthetic clusters will be constructed with different initial conditions implying different distribution for the inferred NFW profiles from dynamics, $(r_{200}^{D},r_{s}^{D})$. The resulting distributions from the dynamics will be in principle different from those from lensing, hence sourcing gravitational slip under equation (3.11). The derived constraints on $\eta$ can then be mapped to specific theories of gravity. Below, we explain the simulations and statistical analysis in detail, before we dwell upon the predicted constraints.

4 Galaxy cluster simulations and statistical analysis

Predicting the constraints on $\eta$ from current or future galaxy cluster observations requires a sample of galaxy clusters along with predictions for the dynamical and lensing mass of each member of the sample according to equation (3.11). The principal aim of our analysis will be to estimate the number of reliable lensing and dynamics mass reconstructions, which will be available from upcoming imaging and spectroscopic surveys, needed in order to constrain $\eta$ at the percent level of statistical accuracy. For this purpose, we will neglect the impact of systematic effects in mass profile determinations, i.e. deviation from dynamical relaxation and spherical symmetry, as well as errors introduced by non-suitable parametrisations of the density profile. Therefore, we will account only for the statistical power of the presented method. An extensive analysis of the impact of systematics in constraining modified gravity models with cluster mass profiles will be performed in a forthcoming work (Pizzuti et al., in prep.).

The analysis carried out in [24] on the X-ray selected CLASH cluster MACS 1206 has shown that accurate dynamics and lensing mass profile reconstructions allow to constrain $\eta(r)$ down to $30\%$ statistical uncertainties with a single galaxy cluster when the systematic effects are under control.

In reconstructing a cluster mass profile from the dynamics of the member galaxies we proceed in two steps as follows: we first generate a synthetic sample of 15 equilibrium configurations of isolated, self-gravitating systems populated by collisionless point-like tracers according to a spherical NFW profile. Each synthetic cluster is constructed with a different characteristic mass and scale, corresponding to different choices for the $r_{200}$ and $r_{s}$ parameters in the NFW profile, covering two orders of magnitude in cluster masses (see Table 4.2). An interesting outcome of our analysis, as we will see in Section 5, is that the constraints do not show an appreciable dependence on halo masses and concentrations. The velocity dispersion of the tracers is assigned following the solution of the Jeans equation, as explained in Section 4.1. The particular number of 15 synthetic halos is chosen as a fair compromise between computing time and a statistically representative sample of cluster masses and concentrations range. The second step, discussed in Sec. 4.2, consists in inputing positions and velocities of the halo’s particles into the MAMPOSSt procedure of Ref. [15], to perform a maximum likelihood fit of the cluster mass profile based on the tracers velocity field. For the lensing part we simulate the probability distribution $P_{L}(r_{s},r_{200})$ relying on the results of the analysis of Ref. [16], which derived the strong+weak lensing mass profile for the galaxy cluster MACS 1206. The profile is parametrized as a NFW model with $\sim 6\%$ and $30\%$ uncertainties on the mass profile parameters $r_{200}$ and $r_{s}$ respectively. Moreover, it was found that the joint probability distribution of the NFW parameters is almost Gaussian with a non-zero covariance between $r_{s}$ and $r_{200}$.

We thus use a bivariate Gaussian distribution,

[TABLE]

[TABLE]

assumed to be centered on the best fit values of the NFW parameters $\bar{r}_{s},\bar{r}_{200}$ given by the dynamics reconstruction, as we will explain in Section 5. Here $\rho$ indicates the correlation. In general, the virial radius $r_{200}$, which is related to the total cluster mass, can be constrained much better than the characteristic radius of the halo mass profile, expressed in terms of $r_{s}$. As shown e.g. in Table 2 of Ref. [17], typical uncertainties on the scale radius given by a lensing probe are of the order of $\sim 30\div 40\%$, while $r_{200}$ can be recovered up to $\sim 5\div 10\%$. For our reference analysis we consider $\sigma_{r_{200}}=0.07\times r_{200}$ and $\sigma_{r_{s}}=0.30\times r_{s}$, while for the correlation in $P_{L}(r_{s},r_{200})$ we adopt the value $\rho=0.67$, found by fitting a bivariate Gaussian on the posterior distribution of Ref. [16]. We further discuss the robustness of our results when varying both correlation and variances in the lensing distribution.

4.1 Synthetic galaxy clusters

Here we discuss the procedure for the generation of the simulated clusters and we list the main properties of the 15 produced synthetic halos we will use in our analysis. Under the assumption of dynamical relaxation, member galaxies in clusters act as collisionless tracers of the total gravitational potential and their dynamics is governed by the Jeans equation. In order to avoid systematic effects induced by lack of equilibrium in the mass profile determination, we generate each halo as an isolated system of collisionless particles assuming spherical symmetry, dynamical relaxation and 3D-Gaussian distribution of the velocities of the tracers (galaxies). Given a $\Lambda$CDM background with $H_{0}=70\,\textrm{km}\,\textrm{Mpc}^{-1}\textrm{s}^{-1}$, $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$444Note that the background cosmology enters only in the definition of $M_{200}$ for each cluster., we populate each halo at redshift $z=0$ according to a NFW density profile, eq. (3.6), up to $\sim 9$ viral radii.

It is natural to expect that the accuracy of the cluster’s reconstructed dynamical mass will depend on the total number of its galaxy members within $r_{200}$. For this purpose, we perform our analysis for two different choices in this regard. In the first case, all 15 synthetic clusters are generated with the same number of tracers equal to 1000, an optimistic, although not unrealistic choice, given the expectation of the number of cluster galaxies with accurate measured spectroscopic redshift from present and future surveys. In the second case we assigned the number of tracers assuming that it is dependent on the cluster mass, requiring a minimum of 300 particles within $r_{200}$ for the less massive halo in the sample. In particular, for the total number of tracers within $r_{200}$ of the $i$-th cluster $N^{(i)}(r_{200})$ we use

[TABLE]

We choose $N_{c}=375$ in order to guarantee a maximum number of tracers less than 1400 for the most massive cluster. It is worth to notice that the choice of eq. (4.2) is totally arbitrary and still an optimistic simplification, since the number of spectroscopically confirmed member galaxies in clusters strongly depends on several factors, such as observational time, completeness of the sample (see e.g. Ref. [22]) and contamination from field galaxies. The resulting $N_{g}(r_{200})$ are displayed in the fifth column of Table 1.

The end result of the sample consists of 15 objects with masses spawning a range from $[5.1\cdot 10^{13}\,M_{\odot}\div 2.9\cdot 10^{15}\,M_{\odot}]$. For each halo, we determine the scale radius $r_{s}$ according to the $c$-$M$ relation (extrapolated at $z=0$) of Ref.[25] 555See also Ref. [26]., which analyzed 19 of the X-ray selected clusters of the CLASH sample at $z\sim 0.19-0.89$ to find

[TABLE]

with

[TABLE]

We have randomly sampled the value of the parameters $A,B$ and $K$ for every halo assuming that they are Gaussian-distributed around the central value and standard deviation of (4.4). The characteristic end-values for each synthetic cluster are summarized in Table 1.

In order to better explore the dependency of the derived constraints on the density of tracers used in the dynamics analysis, we further consider an additional synthetic cluster with $r_{200}=1.96\,\mbox{\rm Mpc}$, $r_{s}=0.27\,\mbox{\rm Mpc}$ generated six times by varying the number of particles $N_{g}(r_{200})$ from $50$ to $2000$.

To each particle at a radial distance $r$ from the center we assign a velocity whose components in spherical coordinates are Gaussian-distributed with a squared dispersion $\vec{\sigma}^{2}_{\mathbf{r}}(r)=\{\sigma^{2}_{r}(r),\sigma^{2}_{\theta}(r),\sigma^{2}_{\phi}(r)\}$, given by the solution of the (spherical) Jeans equation, eq. (3.1). Note that the velocity anisotropy profile $\beta(r)$ is in principle unknown, and one of the major sources of systematics in the reconstruction of cluster mass profile with the analysis of the dynamics of member galaxies. Different methods have been developed in the literature to reconstruct $\beta(r)$, e.g through fitting suitable parametrised profiles, or non-parametric techniques based on inversions of the Jeans equation (e.g. refs. [27, 28]), however, the latter generally requires additional information and assumptions. In our case, we will use a fixed parametric form for $\beta(r)$, namely the Tiret model ”T” of Ref. [29],

[TABLE]

where the normalization of $\beta_{\infty}$ identifies the value of the anisotropy at large radii, while $r_{c}$ is a characteristic scale radius. The Tiret profile is a generalized version of the Mamon and Lokas profile [30], which has been shown to provide a good fit to the average cluster-size halos anisotropy profile over a set of cosmological simulations (see e.g Ref. [31]). Here, we will assume that $r_{c}$ coincides with the scale radius of the mass profile $r_{s}$, and we will also fix $\beta_{\infty}=0.5$ to generate our sample. With the above, the solution for the radial dispersion $\sigma^{2}_{r}$ is given by [30]

[TABLE]

which is connected to the angular one through

[TABLE]

Since we are assuming spherical symmetry we impose $\sigma^{2}_{\phi}\equiv\sigma^{2}_{\theta}$.

4.2 Reconstructing the clusters dynamical mass: the MAMPOSSt method

Armed with the previously constructed synthetic clusters, each with a given $(r_{s},r_{200})$ and a velocity distribution for its tracers, our next step is the reconstruction of what would be the cluster’s mass profile from the kinematical data.

To this purpose, we employ the MAMPOSSt (Modeling Anisotropy and Mass Profiles of Observed Spherical Systems) procedure developed by Ref. [15], which provides a numerical framework for deriving the inferred cluster mass profiles from the analysis of the dynamics of the member galaxies. The code performs a maximum likelihood-fit to the distribution of the positions and velocity field of the galaxies in projected phase space $(R,v_{z})$, where $R$ is the projected radius of each tracer from the cluster center and $v_{z}$ the velocity along the line of sight (los), in the cluster rest-frame. As an example, in the left panel of Figure 1 we show the projected phase space of the $12^{th}$ synthetic cluster in the sample, generated with $1024$ particles within $r_{200}$. The radius $R=r_{200}$ is identified by the red vertical line, while the blue points indicate all the cluster members lying in projection within $r_{200}$.

Given a parametric form of the mass profile $M(r)$ and of the velocity anisotropy profile $\beta(r)=1-{\sigma^{2}_{\theta}}/{\sigma^{2}_{r}}$, the code solves the spherical Jeans equation (3.1) to derive the probability density of finding an object at position $(R,v_{z})$ in the projected phase space as

[TABLE]

Here $N_{p}(R)$ stands for the predicted cumulative number of objects at projected radius $R$ ,and $g(R,v_{z})$ is the surface probability density of observed objects, which in the case of a 3D-Gaussian velocity distribution takes the form

[TABLE]

Given eq. (4.8), the likelihood is simply given by the sum over all the tracers in the projected phase space,

[TABLE]

where as usual, $\vec{\theta}$ stands for the model-parameters vector. The version of MAMPOSSt code used within this paper additionally requires a parametric model of the projected number density profile of the tracers with a characteristic scale radius $r_{\nu}$. Since we are working with collisionless particles, by construction the number density profile $\nu(r)$ scales exactly as the NFW mass profile $\rho(r)$ (i.e. $r_{\nu}\equiv r_{s}$), however, the projected number density profile, needed to compute the probability of eq. (4.8), is obtained in MAMPOSSt by integrating the the 3-dimensional profile along the line of sight, assuming that it extends to infinity. This leads to a value of $r_{\nu}$ which can be slightly different from $r_{s}$ and for this reason we fit the projected number density profile from the phase space of each synthetic cluster and use the best fit values of $r_{\nu}$ as the input for the MAMPOSSt analysis (see the fourth column of Table 1). This is also more reliable from an observational point of view: in a real cluster it is not guaranteed that the distribution of the tracers (i.e. the member galaxies) scales in the same as the total density profile (see e.g. Refs. [32, 33]).

5 Results

We now proceed with the results of our analysis aimed at constraining the gravitational slip $\eta$ from the sample of simulated synthetic clusters described in Section 4. In particular, we will apply the MAMPOSSt procedure to each halo in order to simultaneously constrain the mass profile parameters $(r_{s},\,r_{200})$ and the velocity anisotropy parameter $\beta_{\infty}$, considering all the particles with a projected radius $R\in[0.05\,\mbox{\rm Mpc},\,r_{200}]$ in the fit. This choice has been made to be consistent with real observations, although the range on which we perform the fit could be arbitrary in this ideal case. Indeed, for an observed cluster we cannot assume the validity of the Jeans equation beyond $r_{200}$, while at very small radii the dynamics is dominated by the presence of the brightest central galaxy (BCG), which is the dominant galaxy that sits at the bottom of the potential well of relaxed galaxy clusters (see e.g. Ref. [22]). The final likelihood distribution of dynamics $\mathcal{L}^{\text{dyn}}=\mathcal{L}^{\text{dyn}}(r_{200}^{\text{dyn}},r_{s}^{\text{dyn}})$ is then obtained by marginalizing over $\beta_{\infty}$. In Figure 2 we show the resulting $1\sigma$ and $2\sigma$ iso-probability contours for all the synthetic clusters obtained for what will be our reference case with $N_{g}(r_{200})=1000$ (left plot), and the case with a variable number of tracers within $r_{200}$ (right plot).

It is worth to point out that the best fit values $\bar{r}_{s},\bar{r}_{200}$ given by the MAMPOSSt analysis are not expected to coincide with the true values listed in Table 1, even if these are always included within the $1\sigma$ region. This fact is not surprising, since each synthetic projected phase space represents only one possible statistical realization of a halo with given mass, concentration and anisotropy (we stress again that velocities of the tracers are randomly generated over a Gaussian distribution with given $\vec{\sigma}^{2}_{\mathbf{r}}(r)$). For the purposes of our analysis, we assume as fiducial model the GR expectation $\eta=1$, which is obtained by imposing that the lensing probability distribution $\mathcal{L}^{\text{len}}=P_{L}(r_{200}^{\text{len}},r_{s}^{\text{len}})$ peaks at the best fit values $\bar{r}_{s},\bar{r}_{200}$ given by the dynamics, through an appropriate rescaling which preserves the relative errors $(\sigma(r_{s})/r_{s})^{\text{lens}}$ and $(\sigma(r_{200})/r_{200})^{\text{lens}}$.

5.1 Case A: Scale–independent $\eta$

The first part of the analysis considers $\eta$ to be constant at all scales. This implies that in eq. (3.11) the scale radius of the NFW profile should be the same as inferred by a lensing or dynamics probe, i.e. $r_{s}^{\text{lens}}=r_{s}^{\text{dyn}}\equiv r_{s}$. Thus, the expression of $\eta$ as a function of the mass profiles $M^{\text{lens}},\,M^{\text{dyn}}$ reduces to a constant ratio where the only relevant parameters are the virial radii $r_{200}^{\text{len}},\,r_{200}^{\text{dyn}}$,

[TABLE]

where $c^{\text{lens/dyn}}=r_{200}^{\text{lens/dyn}}/r_{s}$ are the lensing and dynamics concentrations respectively. Following the procedure of Ref. [24], each of the $15$ pairs of joint dynamics-lensing likelihoods are Monte-Carlo sampled for $10000$ sets of values $(r_{s},\,r_{200}^{\text{dyn}}\,,r_{200}^{\text{lens}})$ to obtain the corresponding distribution for $\eta$ according to equation (5.1). The analysis is iterated for the case of synthetic clusters produced with a fixed and a variable number of tracers respectively ($15$ clusters for each case), as explained in Section 4.1.

The results are shown in Figures 3 and 4 . In particular, Figure 3 shows the $1\sigma$ and $2\sigma$ constraints on $\eta$ from each individual cluster, while Figure 4 shows the combined constraints from many clusters. For $N>15$ we extrapolate the predictions on $\eta$ by naively multiplying several times the 15 cluster-combined distribution. An important feature shown in Figure 3 is that the constraints on $\eta$ are practically insensitive to the choice of the member galaxies used within the MAMPOSSt fit, masses and concentration of the cluster, in agreement with the discussion of Ref. [15]. From an observational viewpoint, this is particularly encouraging, since it broadens the range of clusters that may be used from current or future surveys. Furthermore, from Figure 4 it can be seen that the error on $\eta$ obtained by combining together the distributions of many clusters ($N_{\rm clusters}$) scales approximately as the expected theoretical behavior $\sim(N_{clusters})^{-1/2}$ (black dashed line in all panels). We predict that, for our reference analysis (top left panel) the combined analysis of 15 clusters lead to a constraint

[TABLE]

which corresponds to an average statistical uncertainty of $5.5\%$. Furthermore, it turns out that with 75 clusters we reach the accuracy of $\sim 2\%$ and $\sim 4\%$ at $1\sigma$ and $2\sigma$ respectively. When varying the number of particles in the dynamics fit (top right panel) the scaling relation mildly deviates from $(N_{\textrm{clusters}})^{-1/2}$, as a consequence of the slightly larger uncertainties on $\eta$ from the analysis of the smallest halos. The effect is however averaged out when increasing $N_{\rm clusters}$ and the combined constraints are identical to the reference case.

In order to investigate the effect on our predictions induced by variations in the lensing information (still relying on the simplified assumption of a Gaussian distribution) we repeat the analysis doubling the uncertainties in $r_{200}$ and changing the value of the correlation $\rho$ in the lensing information. We have considered three alternatives for $\rho$: a totally uncorrelated distribution ($\rho=0$), a small correlation $\rho=0.15$ and $\rho=0.99$. The results are shown in the lower panels of Figure 4 (for the correlation we plot only the cae $\rho=0.15$) . For completeness, it is worth noticing that we have also considered a change in the relative uncertainties of the scale radius $\sigma(r_{s})/r_{s}$ from $0.3\to 0.45$, but the variation has only marginal effects on the final distributions. The marginalized constraints on $r_{s}$ obtained by the MAMPOSSt fit on our synthetic sample are on average of the same order of the one used in the lensing distribution for the reference analysis $\sigma(r_{s})/r_{s}\sim 0.3$. Therefore, modifying the lensing information on $r_{s}$ in practice does not affect the final predictions for $\eta$. We nevertheless emphasize that mild effects can be enhanced or suppressed depending on the adopted value of the correlation $\rho$. All the constraints obtained in this case are listed in the second column of Table 2.

5.2 Case B: Scale–dependent $\eta$

We now consider $\eta\equiv\eta(r)$ to be scale–dependent, i.e. the scale radius of the lensing mass profile ($r_{s}^{L}$), and that of the dynamics mass profile ($r_{s}^{D}$) are independent parameters for each cluster. As shown in Ref. [24], the constraints on $\eta$ depend on the distance from the cluster center, with the scaling of the uncertainties with $r$ determined by the actual values of $r_{s}^{L},\,r_{s}^{D}$. In order to combine the distributions of $\eta$ obtained for each single cluster we considered 100 equally-spaced radial bins in the range $[0.1,1.5]\,\mbox{\rm Mpc}$. The choice of the lower bound is dictated by consistency with real observations – indeed, as already explained earlier, the assumption of the Jeans analysis cannot be met at very small radii where the dynamics is expected to be dominated by dissipative effects. As for the upper bounds, for an observed cluster the analysis cannot generally be performed beyond the virial radius $r>r_{200}$, where lack of dynamical equilibrium as well as contamination from the large scale structure introduces systematic effects in both lensing and dynamics mass profile reconstructions. We considered the value $r_{max}=1.5\,\mbox{\rm Mpc}$ as a fair compromise that accounts for the different scales of the mass profiles. In each radial bin we Monte-Carlo sampled the joint probability distribution $(r_{s}^{L},r_{200}^{L},rs_{D},r_{200}^{D})$ for 10000 trials, and combined the resulting likelihoods for $\eta(r)$ to get the final probability distribution $P(\eta,r)$. As before, for $N_{cluster}>15$ we extrapolated the behaviour of the constraints by simply multiplying different times $P(\eta,r)$ by itself.

In Figures 5, 6 we show the results for the reference analysis of $N_{g}=1000$, and for the case of variable number of tracers in the MAMPOSSt fit respectively. In each plot the left panel refers to one (inner red shaded area) and two (outer red shaded region) $\sigma$ bands for the 15 cluster-combined distributions while the right panel shows the scaling of the $2\sigma$-errors when increasing $N_{cluster}$ from 15 (red, outer) to 75 (gray, inner). The constraints for all the sub-cases analysed in the scale-dependent configuration, evaluated at the radius $r=1.5\,\mbox{\rm Mpc}$, are summarized in the third column of Table 2. As now $r_{s}^{D}\neq r_{s}^{L}$, we consider also the variation of the uncertainties in the lensing scale radius $r_{s}^{L}$ from 30% to 45% which produces an observable effect on the bounds on $\eta$.

5.3 The effect of the density of member galaxies

Realistic surveys cover a wide range of clusters, each with a different number of tracers depending on the details of its evolution and environment. In this regard, an important question regards the dependence of constraints on the density of tracers within the cluster. Here, we will provide an answer considering the case of a scale-independent $\eta$, and we will show that the effect is practically negligible above a given density of tracers. For this purpose, we generated 6 realisations of a new cluster, characterised by fixed $r_{200}=1.96\,\mbox{\rm Mpc}$, $r_{s}=0.27\,\mbox{\rm Mpc}$, and a Tiret velocity anisotropy profile with $\beta_{\infty}=0.5$, increasing the number of particles enclosed within $r_{200}$ in each realisation, considering the cases $N_{g}(r_{200})=50,100,200,500,1000,2000$. Notice that we chose to keep the lensing probability distribution identical to the reference case.

In Figure 7 we show the $1\sigma$ and $2\sigma$ error bars (left panel), along with the symmetrised $1\sigma$ and $2\sigma$ errors (right panel) of a scale–independent $\eta$ as a function of $N_{g}(r_{200})$. It appears that above $\sim 100$ tracers within $r_{200}$, the goodness of the resulting bounds on $\eta$ increases only marginally, passing from $27\%$ uncertainties at $1\sigma$ for $N(r_{200})=100$, to $21\%$ for $N(r_{200})=2000$; the effect is almost irrelevant for $N(r_{200})>500$.

Extrapolating the result of $\eta(N_{g}=100)$ with the naive theoretically expected scaling $\sim N_{cluster}^{-1/2}$, we obtain

[TABLE]

This shows that very stringent constraints on the gravitational slip can be obtained with a reasonably small number of clusters for which a hundred measured spectroscopic redshifts are available, together with reliable lensing mass profile reconstructions. It is worth pointing out that the analysis carried out in this paper relies on the assumption that the best fit values of $r_{s},\,r_{200}$ from the simulated lensing distribution always coincides with the MAMPOSSt best fit in order to guarantee that the resulting distribution of $\eta$ is centered over the GR expectation $\eta=1$. However, one should take into account that, even in this ideal situation where all the systematics are neglected, the smaller is the number of tracers in the MAMPOSSt analysis, the larger is the difference between the true value of the mass profile parameters and the one fitted by MAMPOSSt (the so called bias parameter, see e.g. Ref. [15]). In Figure 8 we plot the error bars on $\eta$ obtained by centering the lensing distribution on the true values $r_{200}=1.96\,\mbox{\rm Mpc},r_{s}=0.27\,\mbox{\rm Mpc}$ instead of on the MAMPOSSt best fit values. As expected, the bias in $\eta$ tends to be larger for smaller number of tracers, even if GR is always included within $1\sigma$.

6 Summary and discussion

The gravitational slip parameter $\eta\equiv\Phi/\Psi$ is of paramount importance for tests of gravity at cosmological and astrophysical scales. While a possible future detection of $\eta\neq 1$ would be the smoking gun for a departure from the standard gravity paradigm of GR at late-times in the cosmological evolution, a non-detection would disfavor large families of modified gravity scenarios such as the popular $f(R)$ scalar-tensor theories. Here, we elaborated on the potential of galaxy cluster observations to constrain $\eta$, through a joint analysis of dynamics and lensing information. As we discussed in detail in Section 4.1, our ground for the analysis of dynamics of cluster galaxies was based on a sample of synthetic (simulated) clusters, whose dynamical mass profile were subsequently derived by means of the MAMPOSSt numerical framework, in turn to be combined with the lensing one using the available information from the analysis of the cluster MACS 1206.

It is important to remind that our reconstruction of $\eta$ did not invoke the assumption of a particular gravity models or parametrisation. In our view, $\eta\neq 1$ is sourced through a difference in the respective values of the NFW parameters $(r_{200},r_{s})$ when inferred from dynamics and lensing respectively. This phenomenology inclined view is general and allows to translate the derived constraints on $\eta$ to constraints on the theory space of specific models. Nevertheless, we inevitably had to rely on assumptions such as virialization, spherical symmetry (see below) and a NFW parametrisation for the matter density. In addition, the analysis was performed at a fixed redshift $z=0$, which prevented us to capture any possible time-dependence of $\eta$. As concerns the lensing distribution, it is worth to notice that we have performed our analysis under the simplified assumption of Gaussianity for the lensing probability distribution $P_{L}(r_{s},r_{200})$ (even if well-motivated by the reference observational data of MACS 1206). The effect of a variation of the lensing information such as the relative errors and correlation was investigated to find that it has a mild, yet sizable effect on the derived constraints as can be seen from Table 2.

We can summarise our main results as follows:

1. For a scale-independent $\eta$, and assuming a fiducial with $\eta=1$, we can reach an accuracy of $5\%$ ($9\%$) at $1\sigma$ ($2\sigma$) with $15$ clusters, which is brought down to $2\%$ and $4\%$ respectively when extrapolated to $75$ clusters. The explicit results are also summarized in Table 2, and Figures 5 and 6.

2. In a similar manner, for the more realistic case of a scale-dependent $\eta=\eta(r)$, when evaluated at the reference scale $r=1.5\,\text{Mpc}$ and accounting for $15$ clusters we find an accuracy of $\leavevmode\nobreak\ 10\%$ at $1\sigma$ and $21\%$ at $2\sigma$, which decreases to $\sim 4\%$ and $8\%$ when extrapolating to $75$ clusters. Of course, it should come as no surprise that the constraints are less optimistic compared to the simpler scale-independent case.

3. Optimisation of future data analyses in this context requires an understanding of the dependence of the constraints on factors such as the clusters masses and the number of galaxy members in the clusters. Our investigation revealed that the cluster mass (parametrised through $r_{200}$ in the NFW profile) has a rather mild effect on $\eta$. This is best seen from the constraints based on the individual clusters as a function of the cluster mass, as it is shown in Figure 3. According to the procedure explained in Section 4.1 and shown in Table 2, varying the number of tracers in the clusters according to the halo total mass has either none or a very mild quantitative effect on the constraints. To better understand this behavior, we have further investigated the effect of increasing the density of tracers used in the dynamics fit for one single cluster. The results of this exercise show that after $N_{g}(r_{200})\sim 10^{2}$ the effect on the constraints for a scale–independent $\eta$ is moderate, becoming irrelevant for $N_{g}(r_{200})\gtrsim 500$. (See Figures 7 and 8.)

This paper has been devoted to the study of the constraining power of galaxy cluster mass profiles as a probe for modified gravity. It is worth stressing out again that the analysis carried out addressed the problem only from a statistical point of view, neglecting the effect of possible systematics. In principle, several issues can introduce biases in the determination of the gravitational slip – for example, deviations from spherical symmetry and dynamical relaxation, which are the main assumptions of our method, can be sources of differences between the mass profile inferred from lensing and dynamics of the member galaxies respectively, giving rise to spurious detections of gravitational slip, $\eta\neq 1$. Moreover, as already mentioned in Section 4.1, the velocity anisotropy profile $\beta(r)$ is generally an unknown in analyses of cluster dynamics. In the parametric approach implemented in our version of the MAMPOSSt code, the choice of the anisotropy model introduces additional systematics that should be taken into account. A similar statement can be made concerning the parametrisation of the mass profile, as discussed in Ref. [24]. An additional source of systematics not considered here is the presence of interlopers in the dynamical analysis, namely galaxies that are projected in the cluster area but are not dynamically bound to the cluster.

The large amount of data for several hundred clusters that will be provided by current and future imaging and spectroscopic surveys, can allow for determinations of the gravitational slip parameter down to the percent level. This naturally calls for an accurate calibration of the impact of systematics and understanding of the relevant statistics. Such constraints are complementary to other cosmological probes such as the Cosmic Microwave Background, gravitational waves and electromagnetic counterpart detections, cluster number counts and redshift space distortions, allowing to test gravity at cosmological scales down to an unprecedented level of accuracy.

Acknowledgments

I.D.S is supported by ESIF and MEYS (Project CoGraDS – CZ.02.1.01/0.0/0.0/15_003/0000437). S.C. acknowledges grants by CNES and CNRS. L.P. acknowledges S. Borgani and B. Sartoris for useful comments and discussions.