# MG-MAMPOSSt: A code to test modifications of gravity with the dynamics of
galaxy clusters

Lorenzo Pizzuti,1 Ippocratis D. Saltas2 and Luca Amendola3
1 Osservatorio Astronomico della Regione Autonoma Valle d’Aosta, Loc. Lignan
39, I-11020, Nus, Italy
2 CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2,
182 21 Praha 8, Czechia
3Institute of Theoretical Physics, Philosophenweg 16, Heidelberg University,
69120, Heidelberg, Germany
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present an upgraded version of _MG-MAMPOSSt_ , a code to constrain modified
gravity models with the dynamics of galaxy cluster’s members, based on the
original _MAMPOSSt_ method, along with an application to forecasts. Given an
input for the velocity field of member galaxies, the code solves the Jeans’
equation to provide a likelihood analysis for the free parameters of the
gravity theory, that is, the mass profile and modified gravity parameters. The
new version implements two distinct types of gravity modifications, namely
general chameleon (including $f(R)$ models) and Beyond Horndeski gravity
(Vainshtein screening), and is further equipped with a Monte-Carlo-Markov-
Chain module for an efficient parameter space exploration. The program is
complemented by the ClusterGEN code, capable of producing mock galaxy clusters
under the assumption of spherical symmetry and dynamical equilibrium. The code
is written in Fortran and we plan to make it publicly available. Its technical
aspects will be presented elsewhere.

We demonstrate the potential of the method in an ideal situation by analysing
a set of synthetic, isolated and spherically-symmetric dark matter halos,
focusing on the statistical degeneracies between model parameters. Assuming
the availability of additional lensing-like information, we forecast the
constraints on the modified gravity parameters for the two classes of models
presented from joint lensing+dynamics analyses in view of upcoming galaxy
cluster surveys.

###### keywords:

Modified Gravity – Galaxy clusters – Cosmology

††pubyear: 2020††pagerange: MG-MAMPOSSt: A code to test modifications of
gravity with the dynamics of galaxy clusters–MG-MAMPOSSt: A code to test
modifications of gravity with the dynamics of galaxy clusters

##  1 Introduction

Galaxy clusters constitute a very useful probe to test gravity theories,
offering a rich laboratory to study the physics of dark matter, dark energy
and gravity. Combinations of spectroscopic and lensing observations allow for
the reconstruction of the cluster’s mass profile which in turn provide
valuable information about the local gravitational potentials. This idea has
formed the basis for a multitude of phenomenological tests for dark energy
theories beyond General Relativity (GR). The most popular theoretical
frameworks in this context have been scalar-tensor theories, which extend GR
through a new dynamical scalar field. Their theory space has seen strong
constraints after the recent gravitational wave observations (e.g. Ezquiaga &
Zumalacárregui (2017); Creminelli & Vernizzi (2017); Sakstein & Jain (2017);
Baker et al. (2017); Dima & Vernizzi (2018); Kobayashi & Hiramatsu (2018);
Amendola et al. (2018)). The scalar field propagates a new force which
modifies the usual Newtonian one at local scales leaving a phenomenological
footprint on the dynamics and lensing of galaxy clusters. The family of
scalar-tensor theories includes the conformally-coupled theories such as
Brans-Dicke or $f(R)$ gravity, or the more general Horndeski, beyond Horndeski
and the recently introduced DHOST theories (Zumalacárregui & García-Bellido,
2013; Ben Achour et al., 2016; Langlois et al., 2017; Langlois, 2018;
Kobayashi, 2019; Amendola et al., 2020). Previous works have placed strong
constraints on $f(R)$ gravity through the abundance of galaxy clusters
(Schmidt et al., 2009; Rapetti et al., 2010; Lombriser et al., 2012; Rapetti
et al., 2011; Ferraro et al., 2011; Cataneo et al., 2015, 2016b), or combining
reconstructions of the cluster’s thermal and lensing mass (Terukina &
Yamamoto, 2012a; Terukina et al., 2014a; Wilcox et al., 2015a; Sakstein et
al., 2016; Salzano et al., 2016, 2017b). What is more, the combination of
lensing and dynamics reconstruction of a galaxy cluster mass profile can
provide a powerful test based on the gravitational slip parameter $\eta$. This
was pursued in Refs. Pizzuti et al. (2016); Pizzuti et al. (2019); Pizzuti et
al. (2017) which forecasted galaxy cluster constraints on the gravitational
slip following a simple phenomenological approach.

Building on previous works (eg. Pizzuti et al. (2017)), here we present the
upgraded _MG-MAMMPOSSt_ , a new version of the _MAMPOSSt_ code of Ref. Mamon
et al. (2013) developed to determine galaxy clusters mass profiles by
analysing the dynamics of member galaxies in the cluster. Given an input of
line of sight velocities and projected positions of member galaxies, coming
either from simulated or real clusters, the code solves the Jeans equation to
reconstruct the local gravitational potential and the velocity anisotropy
profile, under the assumptions of dynamical relaxation and spherical symmetry.
A key output is a likelihood analysis for the free model parameters in the
mass profile function from dynamics. In a modified gravity scenario, the
gravitational potential receives a contribution from the additional degrees of
freedom, producing an effective mass profile which differs from GR. The first
version of _MG-MAMPOSSt_ has been presented in Pizzuti et al. (2017) for the
case of linear $f(R)$ gravity, and applied to data of two galaxy clusters
analyzed within the CLASH-CLASH/VLT collaborations. In the current version, we
have extended the code’s capability by implementing two popular and viable
modifications of gravity, chameleon and Beyond Horndeski gravity. These models
represent two distinct screening mechanisms, namely the chameleon and
Vainshtein one. They are characterized by their own phenomenological impact on
the inferred mass profile from kinematical data and require a specific study
of their behaviour in this context.

The _MG-MAMPOSSt_ method is complemented by ClusterGEN, a Python3 code capable
of generating synthetic galaxy clusters in a configuration of dynamical
equilibrium, which can be used to perform tests of _MG-MAMPOSSt_ and
forecasting analyses neglecting all the systematic effects.

Our main goal is to showcase the code for the aforementioned class of modified
gravity models, and present an analysis of statistical degeneracy between
model free parameters, which will be essential when testing against real data.
Therefore, we will focus on a catalog of synthetic dark matter halos produced
by means of ClusterGEN. By assuming the availability of additional information
on the galaxy cluster mass profiles, such as that provided by reliable
strong+weak lensing analyses, we also forecast future constraints on the
modified gravity parameters, through a proper implementation of the modified
gravity dynamics in the Jeans equation. In a future paper, we will present an
analysis of real galaxy cluster data.

The paper is structured as follows: In Section 2 we present a brief summary
about the two modified gravity models of this work. In Section 3 we present
the _MG-MAMPOSSt_ code, its features and the new modules implemented with
respect to the original version of Ref. Mamon et al. (2013). We further
introduce the main properties of ClusterGEN, the generator of synthetic halos
accompanying _MG-MAMPOSSt_. Sections 4 and 5 are dedicated to the applications
of _MG-MAMPOSSt_ over a set of simulated Dark Matter halos in the case of
Beyond Horndeski and chameleon gravity respectively. In Section 6 we discuss
our results and main conclusions.

##  2 Vainshtein and chameleon screened theories

The vast majority of scalar-tensor theories rely on a screening mechanism to
recover agreement with local gravity tests at short scales. The dominant
screening mechanisms are the so–called chameleon and Vainshtein mechanism (see
e.g. Burrage & Sakstein (2016); Koyama & Sakstein (2015)), and the current
version of the _MG-MAMPOSSt_ code implements families of theories associated
with both mechanisms.

In the presence of a screening mechanism, the fifth force is allowed to
operate at sufficiently large scales, but gets suppressed locally. Both
chameleon and Vainshtein mechanisms rely on the non-linear interactions of the
scalar field to weaken the coupling between the scalar field and matter.
However, each of them operates differently – whereas chameleon screening needs
the non-linear potential interactions of the scalar, Vainshtein relies on
higher-order derivative self-interactions of the field. Thankfully, either one
does not completely eliminate the trace of the fifth force, but leaves a
characteristic imprint at local scales. Typical examples of theories
exhibiting chameleon screening are $f(R)$ models or scalar-tensor models of
the Brans-Dicke type. On the other hand, Vainshtein screening is naturally a
feature of the more general interactions in Horndeski, Beyond Horndeski and
DHOST scalar-tensor theories.

The different phenomenology of chameleons and Vainshtein is due to the
structurally distinct way the resulting fifth force depends on the source’s
density or its gradients. In theories with chameleon screening the fifth force
is sourced by the scalar potential $\phi$,

$\nabla^{2}\phi=\frac{\partial
V(\phi)}{\partial\phi}+\frac{\mathcal{Q}\rho_{m}}{M_{p}}e^{\mathcal{Q}\phi/M_{p}},$
(1)

with $\mathcal{Q}$ a dimensionless coupling constant111In the literature the
coupling constant is often indicated by $\beta$. Throughout this paper we use
$\mathcal{Q}$ instead to avoid confusion with the anisotropy profile parameter
(see Section 3)., $\rho_{m}$ the density of the matter field and $V(\phi)$ is
a monotonic function of the scalar field. A typical form of the potential
assumes a power-law $V(\phi)=\Lambda^{4+n}\phi^{-n}$ with $n,\Lambda$
constants. The above equation captures a wide range of theories depending on
the value of the coupling $\mathcal{Q}$. In particular,
$\mathcal{Q}=1/\sqrt{6}$ for $f(R)$ gravity. Deep within the massive source,
we have $\nabla^{2}\phi\approx 0$, making the scalar field to reach a minimum,
$\phi_{0}$. Therefore, the fifth-force is screened. Towards the outskirts of
the source, the gradient of the field grows, and leads to a fifth-force given
by

$\displaystyle F_{\phi}=-\frac{\mathcal{Q}}{M_{p}}\frac{d\phi}{dr}.$ (2)

For $f(R)$ gravity the additional degree of freedom is expressed in terms of
the scalaron field
$|f_{R}|\equiv\frac{df}{dR}=\exp(-\sqrt{\frac{2}{3}}\frac{\phi}{M_{p}})-1.$.
Equation (1) requires the input of a profile for the matter density $\rho$. In
this work, we will model the halo mass profiles assuming a NFW profile
(Navarro et al., 1997) as

$\rho_{m}(r)=\frac{\rho_{s}}{r/r_{s}(1+r/r_{s})^{2}},$ (3)

with $\rho_{s}$ the typical density and $r_{s}$ the scale radius at which the
logarithmic derivative of the profile equals $-2$. For this model, following
e.g. Ref. Terukina & Yamamoto (2012b), it is possible to derive an analytical
approximation for the radial profile of $\phi(r)$, which matches the two
limiting solutions inside (where $\nabla^{2}\phi\approx 0$) and outside the
object (where $\partial V(\phi)/\partial\phi$ is subdominant),

$\phi(r)=\begin{cases}\phi_{s}[r/r_{s}(1+r/r_{s})^{2}]\equiv\phi_{int}(r)&r<S\\\
-\frac{\mathcal{Q}\rho_{s}r_{s}^{2}}{M_{p}}\frac{\ln(1+r/r_{s})}{r/r_{s}}-\frac{C}{r/r_{s}}+\phi_{\infty}\equiv\phi_{ext}(r)&r>S.\end{cases}\
$ (4)

In the above equation
$\phi_{s}=(n\Lambda^{n+4}M_{p}/\mathcal{Q}\rho_{s})^{1/(n+1)}$ is the value of
the field in the interior region, $\phi_{\infty}$ the background value and $C$
an integration constant. The _screening radius_ $S$ represents the transition
scale between the two regimes. $C$ and $S$ can be determined requiring that
$\phi_{int}$ and $\phi_{ext}$ and their derivatives match at $S$. In the
halo’s inner core the field is strongly suppressed with respect to its ambient
value $\phi_{s}\ll\phi_{\infty}$, and we can approximate the interior solution
to be negligible $\phi_{int}\simeq 0$. In this case, for $C$ and $S$,

$S=\frac{\mathcal{Q}\rho_{s}r_{s}^{3}}{M_{p}\phi_{\infty}}-1,\;\;\;C=-\frac{\mathcal{Q}\rho_{s}r_{s}^{2}}{M_{p}}\ln(1+S/r_{s})+\phi_{\infty}S/r_{s}.$
(5)

For beyond Horndeski and DHOST theories which agree with gravitational wave
observations, it has been shown that the Vainshtein mechanism is partly broken
within massive sources. This leads to a different fifth force term, which does
not depend only on the source’s density but also on its gradient (Kobayashi et
al., 2015; Crisostomi & Koyama, 2018; Dima & Vernizzi, 2018),

$\frac{d\Phi(r)}{dr}=\frac{GM(r)}{r^{2}}+\frac{Y_{1}G}{4}\cdot\frac{d^{2}M(r)}{dr^{2}}.$
(6)

Here, $M(r)$ is the mass of the object enclosed within radius $r$, and $Y$ the
dimensionless coupling of the new force. The gradient of the mass profile is
given as standard by $dM(r)/dr=4\pi\rho_{m}(r)r^{2}$. In a galaxy cluster, the
density $\rho_{m}(r)$ is assumed to be the density of dark matter as before. A
$Y_{1}>0$ ($Y_{1}<0$) weakens (strengthens) gravity, remembering that
$d\rho_{m}/dr<0$ in the interior. Standard gravity is recovered outside the
object. Currently, the upper and lower bound from astrophysical probes come
are of the order $Y<10^{-2}-10^{-1}$ respectively ( Saltas & Lopes (2019);
Crisostomi et al. (2019); Saltas et al. (2018); Sakstein (2018); Ishak (2019);
Creminelli et al. (2018); Babichev & Lehébel (2018); Dima & Vernizzi (2018);
Sakstein et al. (2016); Babichev et al. (2016); Sakstein (2015a, b); Jain et
al. (2016); Sakstein (2013)).

A fundamental difference between the two previous theories is that, whereas
conformally-coupled theories such as $f(R)$ leave lensing unaffected, this is
not any more true within beyond Horndeski and DHOST theories. In the latter
case, the relativistic potential $\Psi$ is also modified explicitly as

$\displaystyle\frac{d\Psi(r)}{dr}=\frac{GM(r)}{r^{2}}-\frac{5Y_{2}G}{4r}\cdot\frac{dM(r)}{dr},$
(7)

with $Y_{2}$ a dimensionless coupling. Lensing is sourced by the combination
$\Phi_{\text{lens}}\sim\frac{1}{2}(\Phi+\Psi)$, and it may provide constraints
on the additional coupling $Y_{2}$ present in these theories (e.g. Sakstein
(2018)). In particular, we can define an effective lensing mass which,
assuming spherical symmetry, is given by

$M_{\text{lens}}=\frac{r^{2}}{2G}\left[\frac{d\Psi}{dr}+\frac{d\Phi}{dr}\right],$
(8)

with $G$ the bare Newton’s coupling. Notice that $M_{\rm lens}$ is sensitive
to both $Y_{1}$ and $Y_{2}$. This will be useful in analyses where dynamical
probes are combined with lensing measurements.

##  3 The _MG-MAMPOSSt_ and ClusterGEN codes

Here we introduce the main features of _MG-MAMPOSSt_ , focusing on the new
modules implemented for modified gravity. We also briefly describe the
ClusterGEN code, which is used to produce the mock clusters for the
applications of _MG-MAMPOSSt_ presented in Sections 4 and 5.

Figure 1: Diagram which illustrates the structure of the _MG-MAMPOSSt_ code.

### A. MG-MAMPOSSt

_MG-MAMPOSSt_ is a new version of the _MAMPOSSt_ code of Ref. Mamon et al.
(2013) for the dynamical analysis of collisionless, cluster-size systems aimed
at constraining modified gravity models. The method has been first applied in
Ref. Pizzuti et al. (2017) for the simple case of linear $f(R)$ gravity. The
upgraded version of the code implements general non-linear chameleon gravity
and Beyond Hordenski theories.

Main concept: The _MAMPOSSt_ (Modelling Anisotropy and Mass Profile of
Spherical Observed Systems) procedure, originally developed in Mamon et al.
(2013), infers the mass profile and orbit anisotropy profile of galaxy
clusters under the assumption of spherical symmetry and dynamical relaxation.
The code uses the projected information on the cluster member galaxies, which
are considered as collisionless tracers of the gravitational potential, i.e.
the projected radius $R$ from the center of the cluster at which a galaxy is
seen by the observer and the line-of-sight (l.o.s) velocity $v_{z}$, computed
in the rest frame of the cluster. Working in this projected phase space
$(R,v_{z})$, _MAMPOSSt_ performs a Maximum Likelihood fit by solving the
(static) spherical Jeans’ equation,

$\frac{d(\nu(r)\sigma_{r}^{2})}{dr}+2\beta(r)\frac{\nu\sigma^{2}_{r}}{r}=-\nu(r)\frac{d\Phi}{dr}.$
(9)

In the above equation $\nu(r)$ corresponds to the number density profile of
tracers, $\sigma^{2}_{r}$ is the velocity dispersion along the radial
direction, $\beta\equiv
1-(\sigma_{\theta}^{2}+\sigma^{2}_{\phi})/2\sigma^{2}_{r}$ the velocity
anisotropy profile and $\Phi$ is the total gravitational potential, which can
provide a definition for the dynamical mass profile as

$\frac{d\Phi}{dr}\equiv\frac{GM_{\text{dyn}}(r)}{r^{2}}.$ (10)

Anisotropy profile: The anisotropy profile $\beta(r)$ is not known in
principle and represents one of the major systematics in the dynamical mass-
profile reconstruction. To infer it, both parametric (e.g. Biviano et al.
(2013); Wojtak et al. (2009)) and non-parametric (e.g. Stark et al. (2019);
Host et al. (2009)) methods has been developed. In _MG-MAMPOSSt_ we adopt a
parametric approach with six different models. However, for illustrative
purposes, we will rely only on the Tiret model of Ref. Tiret et al. (2007),

$\beta_{T}(r)=\beta_{\infty}\frac{r}{1+r_{c}}.$ (11)

$\beta_{\infty}$ is the anisotropy at large radii and $r_{c}$ a scale
parameter which in _MG-MAMPOSSt_ has been set to be equal to the scale radius
of the mass profile $r_{s}$. The Tiret model is a general version of the Mamon
& Lokas profile of Mamon & Łokas (2005), which has been shown to provide an
average good fit to the anisotropy of cluster-size halos from cosmological
simulations (e.g Mamon et al. (2010)). It is worth to point out that the
constraints on MG models exhibit a negligible dependence on the anisotropy
profile chosen in the ClusterGEN generator.

Velocity and number density profiles: In our current version of the code we
assume that the three-dimensional distribution of the velocity field is
Gaussian. In this case, given the expression for the radial velocity
dispersion from Eq. (9) (e.g. Mamon & Łokas (2005)),

$\sigma^{2}_{r}=\frac{1}{\nu(r)}\int_{r}^{\infty}{\exp\left[2\int_{r}^{s}{\frac{\beta(t)}{t}dt}\right]\nu(s)\frac{d\Phi}{ds}ds},$
(12)

and for the velocities along the l.o.s at a given $R$ we have that

$\sigma^{2}_{z}=\left[1-\beta(r)\left(\frac{R}{r}\right)^{2}\right].$ (13)

Note that there is no _a priori_ restriction on the choice of the 3D velocity
profile. However, the tests performed by Mamon et al. (2013) over a set of
halos from cosmological simulation showed that the Gaussian model works
sufficiently well even when the l.o.s velocity distribution exhibits
significant deviations from Gaussianity. The number density $\nu(r)$ can in
general be measured directly from the phase space. As such, it can be excluded
from the fit. We parametrize $\nu(r)$ adopting a projected NFW profile, which
comes out by integrating Eq. (3) along the l.o.s. Since the normalization of
the profile is factored out in Eq. (12) the model is entirely parametrised by
the scale radius $r_{\nu}$. It is worth noticing that in general, the
distribution of the cluster’s members differs from the distribution of the
total matter density (see e.g. Biviano & Salucci (2006); Budzynski et al.
(2012)). Therefore, the scale radii $r_{\nu}$ and $r_{s}$ (of the total mass
profile) need not coincide. For additional information about the anisotropy
models and mass profiles implemented in the original version, see Mamon et al.
(2013).

Code’s output and sampling features: _MG-MAMPOSSt_ constrains the free input
parameters related to e.g the mass profile and/or the gravity theory under
study. The output of the code is the logarithm of the total likelihood as

$-\ln\mathcal{L}_{\text{dyn}}=-\sum_{i=1}^{N}q(R_{i},v_{z(i)}|\bf{\theta}).$
(14)

$N$ is the total number of member galaxies considered in the fit,
$q(R_{i},v_{z(i)}|\bf{\theta})$ is the probability of observing a galaxy with
projected position $R_{i}$ and l.o.s velocity $v_{z(i)}$ given the model(s)
described by the set of parameters $\bf\theta$. The expression of $q$ is
determined by the l.o.s velocity distribution at a given $R$ and the projected
number density profile of the tracers.

In the first version of _MG-MAMPOSSt_ , likelihoods were computed over a
multi-dimensional grid of points in the parameter space. However, for a large
number of parameters ($\gtrsim 4$) the grid is computationally too expensive
when the code runs in Modified Gravity mode. In the new version, we have
implemented a Monte Carlo-Markov Chain (MCMC) module, in which the chain is
based on a simple Metropolis-Hastings algorithm with a fixed-step Gaussian
random walk. The code takes few hours to produce a full output of $\sim
10^{5}$ points for the models currently implemented. The MCMC module can read
multiple data files at the same time and run the chain over the combined
likelihood, a tool which is particularly useful for studying simulated
catalogs to perform forecast analyses.
In Fig. 1 a scheme which summarize the structure of _MG-MAMPOSSt_ is shown.

Gravitational potential: Compared to the original _MAMPOSSt_ code, _MG-
MAMPOSSt_ introduces new parametrisations for the total gravitational
potential in Eq. (12). All the models are based on a NFW mass density profile,
but a work-in-progress next version of the code will extend this to other mass
profiles. The particular gravitational potentials implemented in the current
version are as follows:

1. Linearised Horndeski Gravity: These models correspond to a scalar-tensor subclass of Horndeski theories where the Newtonian potential $\Phi$ can be written in terms of three parameters, namely $h_{1}$, $\mathcal{O}$ and $m$ (see e.g. Amendola et al. (2013)). These parameters are functions of redshift, but they can be considered constant in first approximation. In this model, there is no screening at all. In order to pass the local gravity constraints, either one assumes a very small coupling, or that the baryons are not sensitive to the fifth force.

In particular, assuming the NFW profile for matter density perturbations we
have

$\frac{d\Phi}{dr}=h_{1}\left[\frac{GM_{\text{NFW}}(r)}{r^{2}}+2\mathcal{O}^{2}\frac{d\Phi_{MG}}{dr}\right],$
(15)

where $M_{\text{NFW}}$ is the standard NFW mass profile and $\Phi_{MG}$ is
given by

$\Phi_{\text{MG}}(r)=\frac{2\pi
G\rho_{s}}{r}r_{s}^{3}\left\\{e^{-m(r_{s}+r)}\left[{\rm Ei}(m\,r_{s})-{\rm
Ei}(m(r_{s}+r))\right]\right.,$ $\left.-e^{m(r_{s}+r)}{\rm
Ei}\left[-m(r_{s}+r)\right]+e^{m(r_{s}-r)}{\rm Ei}(-m\,r_{s})\right\\}.$ (16)

The above parametrization has been first introduced in _MAMPOSSt_ and tested
by Ref. Pizzuti et al. (2017). For the particular case of $f(R)$, it is
$2\mathcal{O}^{2}=1/3$, $h_{1}\sim 1$ and $m=const\equiv
m_{f_{R}}=[(f_{R}/f_{RR}-R)/3]^{1/2}$, the mass of the additional scalaron
field, where $f_{R}=df/dR$ and $f_{RR}=d^{2}f/dR^{2}$. Note that $m_{f_{R}}$
is the only free parameter to be constrained.

2. Beyond Horndeski Gravity: The form of the gravitational potential is according to (6). The coupling constant $Y_{1}$ is the additional free parameter, along with the mass profile and velocity anisotropy parameters. With a NFW mass density profile, (6) can be rewritten as

$\frac{d\Phi}{dr}=\frac{GM_{200}}{r^{2}(\log(1+c)-c/(1+c))}\left[\log\left(1+\frac{r}{r_{s}}\right)\right.$
$\left.-\frac{r/r_{s}}{\left(1+r/r_{s}\right)}+\frac{Y_{1}}{4}\frac{r^{2}(r_{s}-r)}{(r_{s}+r)^{3}}\right],$
(17)

where $M_{200}=4\pi r_{s}^{3}\rho_{s}[\log(1+c)-c/(1+c)]$ is the total mass of
a sphere of radius $r_{200}$ enclosing an average matter density 200 times the
critical density of the Universe at that redshift $\rho_{c}(z)$. The parameter
$c=r_{200}/r_{s}$ is the concentration. Eq. (17) can be rewritten in a Poisson
GR-like form by defining an effective dynamical mass such that

$\frac{d\Phi}{dr}\equiv\frac{GM_{\text{dyn}}}{r^{2}}=\frac{G}{r^{2}}(M_{NFW}+M_{1}),$
(18)

where

$M_{\text{NFW}}=M_{200}\frac{\log(1+r/r_{s})-(r/r_{s})(1+r/r_{s})}{\log(1+c)-c/(1+c)}$
(19)

is the standard NFW mass profile in GR and

$M_{1}=M_{200}\frac{Y_{1}}{4}\frac{r^{2}(r_{s}-r)}{(r_{s}+r)^{3}}\times[\log(1+c)-c/(1+c)]^{-1}$
(20)

is the contribution induced by the modification of gravity.

In the top left plot of Fig. 2 we show the relative ratios
$(M_{\text{dyn}}-M_{\text{NFW}})/M_{\text{NFW}}$ of the effective mass
profiles $M_{\text{dyn}}$ of a NFW model in Beyond Horndeski gravity,
characterized by $r_{200}=2.0\,\text{Mpc},\,\,r_{s}=0.3\,\text{Mpc}$ for two
different coupling values, $Y_{1}=0.5$ (blue line) and $Y_{1}=-0.5$ (red
line). Compared to GR, a positive $Y_{1}$ enhances the effect of gravity for
radii smaller than $r_{s}$ while it decreases the gravitational potential at
large $r$. The corresponding relative ratios of the velocity dispersion
profiles $\sigma_{r}(r)$, obtained in _MG-MAMPOSSt_ from Eq. (12) assuming a
Tiret anisotropy model with $r_{\beta}=r_{s}$ and $\beta_{\infty}=0.5$222As
for the number density profile, we still assume a NFW with $r_{\nu}\equiv
r_{s}$, are illustrated in the top right panel of 2.

Figure 2:  Top: Relative ratios of the effective mass profiles (left) and of
the velocity dispersion profiles (right), as a function of the distance from
the cluster center, for a NFW halo with
$r_{200}=2.0\,\text{Mpc},\,\,r_{s}=0.3\,\text{Mpc}$ and with a velocity
anisotropy profile described by a Tiret model with $r_{\beta}=r_{s}$ and
$\beta_{\infty}=0.5$. The profiles are obtained in Beyond Horndeski gravity
with $Y_{1}=0.5$ (red curves) and $Y_{1}=-0.5$ (blue curves). The vertical
brown dashed line in the left plot represents the radius $r_{200}$.
$M_{\text{NFW}}(r)$ is the NFW mass profile in GR. For $Y_{1}>0$ the effect of
gravity is enhanced with respect to GR for radii smaller than $r_{s}$ and at
large $r$ and exactly the opposite happens for $Y_{1}<0$. Bottom: relative
ratios of the effective mass profiles in chameleon gravity with
$\mathcal{Q}=1/\sqrt{6}$, $\phi_{\infty}=10^{-4}$ (red curves),
$\mathcal{Q}=1/\sqrt{6},\,\,2\times 10^{-5}$ (blue curves) and
$\mathcal{Q}=1$, $\phi_{\infty}=10^{-4}$ (gray curves). The left plot shows
profiles derived from a NFW halo with
$r_{200}=2.0\,\text{Mpc},\,\,r_{s}=0.3\,\text{Mpc}$ while the right plot is
for halos with $r_{200}=1.5\,\text{Mpc},\,\,r_{s}=0.5\,\text{Mpc}$. In both
plots the vertical brown dashed line represent the radius $r_{200}$.

3. General Chameleon Screening: For this class of models the gravitational potential is given by

$\frac{d\Phi}{dr}=\frac{GM_{\text{NFW}}(r)}{r^{2}}+\frac{\mathcal{Q}}{M_{p}}\frac{d\phi}{dr},$
(21)

where the chameleon field $\phi$ is given by Eq. (4) under the approximations
described in Section 2. The additional degrees of freedom are now the coupling
constant $\mathcal{Q}$ and the ambient field value $\phi_{\infty}$, which sets
a boundary condition for the field outside the source. As mentioned above,
chameleon $f(R)$ gravity is included as the sub-class obtained for
$\mathcal{Q}=1/\sqrt{6}$.

While for linear Horndeski and Beyond Horndeski gravity the modified gravity
parameters are independent of the size of the halo, in chameleon screening the
suppression of the fifth force is related to the matter density. In
particular, as we will show in Section 5, clusters with smaller mass and
larger $r_{s}$ have a screening mechanism less efficient with respect to more
massive halos.

In the bottom plots of Fig. 2 we show the relative ratio of mass profiles in
chameleon gravity with respect to GR obtained for different values of
$\mathcal{Q}$ and $\phi_{\infty}$. In both plots, solid red and blue lines
refer to an $f(R)$-type model ($\mathcal{Q}=1/\sqrt{6}$) with
$\phi_{\infty}=10^{-4},\,\,2\times 10^{-5}$ respectively. For a large
background field the cluster is not screened and the effective mass profile is
enhanced of a factor $2\mathcal{Q}^{2}$ compared to the standard NFW model.
Decreasing $\phi_{\infty}$ the screening mechanism becomes active and
$M_{\text{dyn}}\equiv M_{\text{NFW}}$ for $r\leq S$. Grey curves shows the
effective mass profiles for a chameleon field $\phi_{\infty}=10^{-4}$ and a
coupling constant $\mathcal{Q}=1$; since the screening effect is connected
also to the coupling with matter, part of the halo is screened while in the
non-screened region the fifth force increases. In the bottom left plot of Fig.
2 the NFW parameters are $r_{200}=2.0\,\text{Mpc},\,\,r_{s}=0.3\,\text{Mpc}$,
while the bottom right plot shows the same chameleon models for a halo with
$r_{200}=1.5\,\text{Mpc},\,\,r_{s}=0.5\,\text{Mpc}$. As expected, the
screening suppression is reduced and the effect of the fifth force is enhanced
also for small field values.

### B. ClusterGEN

The ClusterGEN code, first presented in Ref. Pizzuti et al. (2019) is a
generator for isolated, self-gravitating systems in dynamical equilibrium,
populated by collisionless point-like tracers according to a spherical model.
The code produces mock-data catalogs which serve as input for _MG-MAMPOSSt_.
The main purpose is to test the constraining power of the _MG-MAMPOSSt_ method
for different modified gravity scenarios in ideal conditions where all the
systematics are under control. The synthetic clusters can be built with
different characteristic masses and scales, while the dynamics is assumed to
be governed by the Jeans’ equation (9) with a Gaussian 3D distribution for the
velocity field.

For a given radial number density profile, ClusterGEN distributes the
particles in spherical shells from the cluster center and assigns to each
tracer a velocity whose components, in spherical coordinates, have a squared
dispersion
$\vec{\sigma}^{2}_{\mathbf{r}}(r)=\\{\sigma^{2}_{r}(r),\sigma^{2}_{\theta}(r),\sigma^{2}_{\phi}(r)\\}$.
$\sigma^{2}_{r}(r)$ is given by the integral solution Eq. (12) and
$\sigma^{2}_{\theta}(r)=\sigma^{2}_{\phi}(r)$ is then obtained from the
velocity anisotropy profile as

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

The current version of ClusterGEN works with four different anisotropy models,
namely constant anisotropy $\beta(r)=\beta_{C}$ (which reduces to the
isotropic case when $\beta_{C}=0$), the Mamon & Lokas profile, the Tiret
profile and the modified Tiret profile (see e.g. Biviano et al. (2013)), which
allows for negative anisotropy in the halo center. As for the gravitational
potential, the radial velocity dispersion of Eq. (12) can be obtained both
assuming GR, and with one of the modified gravity models implemented in the
_MG-MAMPOSSt_ code. For the purposes of this paper, we study the constraints
on the modified gravity parameters obtained from _MG-MAMPOSSt_ when the
procedure is applied only on GR mock catalogs. We will leave the analysis of
non-GR synthetic halos for future applications.

The projected phase spaces $(R,\,v_{z})$ which enter the _MG-MAMPOSSt_ code
are then extracted along the chosen projection axis following Eq. (13). In
Fig. 3 we plot two phase spaces obtained for one cluster in our catalog (see
next Section for halo’s parameters details) considering $N=600$ particles
(left) and $N=100$ particles (right) within the radius $r_{200}$.

Figure 3: Projected phase spaces for a standard NFW halo with
$r_{200}=2.0\,\text{Mpc},\,\,r_{s}=0.3\,\text{Mpc}$, assuming a Tiret model
for the velocity anisotropy profile with $r_{\beta}=r_{s}$ and
$\beta_{\infty}=0.5$. Left: 600 tracers considered within $r_{200}$ (vertical
brown dashed line). Right: 100 tracers within $r_{200}$.

##  4 Application I: Beyond Horndeski gravity

In this Section we show the application of _MG-MAMPOSSt_ for the choice of
Beyond Horndeski, to spherically-symmetric and virialised synthetic clusters
generated with the ClusterGEN code as outlined in Section 3. We investigate
the possible constraints on the additional degree of freedom $Y_{1}$ from
dynamics analysis alone and simulating the availability of additional
information such as the probability distribution provided by a strong+weak
lensing probe.

###  4.1 Setup

For Beyond Hordenski models (BH), modifications of gravity are not related to
the size of the halo considered. Thus, all synthetic phase spaces are
generated as different realization of the same massive NFW halo with
$r_{200}=2.0\,\text{Mpc}$, $r_{s}=0.3\,\text{Mpc}$, and with a Tiret
anisotropy profile $\beta(r)=\beta_{\infty}r/(r+r_{\beta})$ with
$r_{\beta}=r_{s}$ and $\beta_{\infty}=0.5$. In view of the proof-of-concept
aim of the paper, this choice has been made to facilitate the combined
analysis of several clusters. Furthermore, in order to avoid spurious
signatures of modified gravity due to statistical fluctuations in the cluster
generation process, each phase space is derived requiring that the best fits
of the standard _MAMPOSSt_ analysis in GR differs from the cluster true values
by a reduced chi square $\hat{\chi}^{2}=0.02$. This way, the obtained
constraints reflect the true statistical degeneracy between the model
parameters.

As a default choice, each synthetic cluster is constructed with 600 particles
(tracers) in the radial range $[0.05\,\mbox{\rm Mpc},r_{200}]$. This choice is
a reasonable, although optimistic, expectation for the number of member
galaxies’ spectroscopic redshifts available for few dozen clusters from
ongoing and next generation surveys. It is, however, an important systematic
factor in realistic observations. We will therefore discuss how our results
are affected when lowering the number of tracers considered in the synthetic
phase spaces. As concerns the projected radial range, we note that for real
data the dynamics of member galaxies in clusters is dominated by the Brightest
Central Galaxy (BCG) in the innermost region, which is thus excluded from the
analysis. Similarly, even for relaxed clusters, departures from dynamical
equilibrium cam become important above $r_{200}$, and the Jeans’ equation is
no more valid. For further discussions see eg. Refs. Biviano et al. (2013);
Pizzuti et al. (2019)

###  4.2 Kinematical analysis

We considered a catalog of 20 synthetic halos generated in GR. This is a fair
expectation of the number of relaxed galaxy clusters for which very high
quality imaging and spectroscopic data could be available from present and
future surveys.

Given the tracers’ phase space $(R,v_{z})$, we first apply the _MG-MAMPOSSt_
procedure to fit for the parameters of the mass profile ($r_{200},\,r_{s}$),
the anisotropy parameter ($\beta_{\infty}$) and the modified gravity parameter
$Y_{1}$. The results for the 2-dimensional marginalized iso-probability
contours and the one-dimensional distributions for each parameters are shown
in Fig. 4 and Fig. 5 respectively, for a single cluster. The green curves in
Fig. 4 represent the iso-contours of $\Delta\chi^{2}=6.1$, while the red inner
lines indicate $\Delta\chi^{2}=2.3$333Note that we choose to compute
confidence intervals using $\chi^{2}$ levels instead of directly integrating
the likelihood. This is because $\chi^{2}$ regions are not sensitive to the
peak of the distribution, but to the overall shape of the likelihood, which in
our case may exhibit irregular features. Moreover, for non-normalizable
distributions, $\chi^{2}$ bounds can be always set, independently of the
integration limits, which is advantageous for the marginalisation process,
although in this case one cannot associate a probability meaning to the
likelihood.. For the velocity anisotropy, from now on we use

$\beta\equiv(1-\beta_{\infty})^{-1/2},$ (23)

as in Ref. Biviano et al. (2013), which is equal to 1 for isotropic
distributions ($\beta_{\infty}=0$). We have marginalized over the logarithm of
the parameters, except for $Y_{1}$ which can assume negative values, in order
to give an equal weight to all values. The blue points as well as the red
lines in Fig. 5 correspond to the true values of the halo’s parameter.

A first important lesson is that the coupling constant $Y_{1}$ is strongly
degenerate with the mass profile parameters – increasing the NFW parameters
$r_{s}$ and $r_{200}$ corresponds to an increase of $Y_{1}$ respectively. This
is expected, since increasing $Y$ weakens the gravitational potential, but so
does increasing the size of the cluster of a given mass. As a consequence, the
effect of a standard Newtonian gravitational potential in a galaxy cluster
phase space can be mimicked to a certain extent by introducing an additional
degree of freedom $Y_{1}$ and suitably adjusting the values of the other
parameters. This aspect needs to be accounted for in the analysis of real
observations when placing constraints on the modified gravity coupling.

Fig. 5 shows the marginalised, 1-dimensional likelihoods for one halo in the
sample. Note that all the distributions of $(Y_{1},\,r_{s},\,r_{200})$ are
shifted towards larger positive values with respect to the expectation, which
is particularly evident for $Y_{1}$. This is a statistical artifact and not a
physical effect – indeed, all the four-dimensional likelihoods exhibit a
maximum in the neighbourhood of the halo’s parameters true values (i.e.
$Y_{1}=0$, $r_{200}=2.0\,\mbox{\rm Mpc}$, $r_{s}=0.3\,\mbox{\rm Mpc}$,
$\beta=1.41$) by construction. However, they also have a second local maximum
at $Y_{1}>1$, due to the strong degeneracy with the mass profile parameters,
which shifts the 1-dimensional distribution when marginalizing. We can always
set an upper limit $Y_{1}<6.02$ at $\Delta\chi^{2}=4.0$ from the marginalised
distribution (corresponding to 95% C.L.). The value has been chosen as the
largest constraint obtained among the 20 marginalized distributions $P(Y_{1})$
of each of the synthetic phase spaces in the sample.

Figure 4: Iso-chi square contours of the marginalised 2-dimensional
distributions for the free parameters in the _MG-MAMPOSSt_ fit, i.e the NFW
profile and modified gravity parameter $Y_{1}$. The results are shown for one
halo in the analysis, but it is representative of the behaviour of all the
cluster sample. Red curves: $\Delta\chi^{2}=2.3$. Green curves:
$\Delta\chi^{2}=6.1$. The blue points indicate the true values of the halo’s
parameters. Figure 5:  Marginalized one-dimensional distributions of the four
free parameters obtained by the _MG-MAMPOSSt_ analysis for the same halo of
Fig. 4. The dashed lines represent the true value of the parameters.

So far we focused on a single cluster. Still working with the synthetic
clusters constructed assuming GR, we combine the individual likelihoods of
each phase space to study the joint distribution as a function of the number
of halos. We write,

$-\log\mathcal{L}_{\text{tot}}=\sum_{i}^{N}\left(-\log\mathcal{L}^{\text{dyn}}_{i}\right),$
(24)

where $N$ is the number of phase spaces in the sum. Remember that we are
considering realizations of halos characterized by the same NFW parameters,
but with distinct projected phase spaces. Thus, a stacked analysis can be done
by means of combining single likelihoods of different (synthetic) clusters as
in Eq. (24). In Fig. 6 we plot five marginalized distributions $P(Y_{1})$,
corresponding to different number of halos $N$ taken into account in Eq. (24).
The curves exhibit an apparent tension with $Y_{1}=0$ for $N>1$. When
increasing the number of halos, the distributions develop a second maximum
near the GR expectation value which, however, remains subdominant with respect
to the first peak. As one may naively expect, this confirms that the
statistical degeneracy, responsible of the spurious maximum in the
marginalized likelihood, cannot be removed just by increasing the number of
halos in the fit. For $N\sim 20$ we may extract $Y_{1}\lesssim 3.9$ at $95\%$
C.L. As for the lower limit, we are marginally close to the the theoretical
lower bound, $Y_{1}>-2/3$ for the combined distribution of $20$ clusters.

The results of our analysis indicate that, in the case of BH gravity, the
constraints obtainable with the _MG-MAMPOSSt_ fit alone are limited by the
statistical degeneracy between model parameters. Other probes, such as X-ray
and lensing measurements, are required in order to remove the unphysical
region of the parameter space associated to large values of
$r_{s},\,r_{200},\,Y_{1}$. Similar conclusions have been already derived in
Ref. Pizzuti et al. (2017) for the application of _MG-MAMPOSSt_ to linear
$f(R)$ gravity. Therefore, although kinematics alone are not able to provide
sufficiently constraining bounds on $Y_{1}$, from the combination of several
data-sets obtained with different methods one may expect to obtain valuable
results. In order to further investigate this point, in the following we
simulate the availability of additional lensing information to understand at
what level of precision the models can be constrained with our method.

Figure 6:  Marginalized distributions $P(Y_{1})$ obtained by combining
likelihoods of single halos according to Eq. (24). Different curves correspond
to different number of phase spaces considered in the combination. The effect
of strong degeneracy between $Y_{1}$ and the NFW parameters $r_{s},\,r_{200}$,
is particularly evident for $N_{halos}\gtrsim 5$, giving rise to a spurious
peak at $Y_{1}>0$. The value of $P(Y_{1})$ corresponding to the second
maximum, near the GR expectation ($Y_{1}=0$), increases with the number of
halos but is always subdominant with respect to the mode of the distribution.

###  4.3 Joint Lensing and dynamics analysis

As already discussed in Ref. Pizzuti et al. (2019), present and upcoming high
precision Strong+Weak lensing measurements can recover the mass profile
parameters $r_{200}$ and $r_{s}$ with uncertainties down to $\sim 7-10\%$ and
$\sim 30\%$ respectively (e.g. Umetsu et al. (2016)). In order to forecast the
constraints obtainable by our method with an additional lensing probe, we
consider a Gaussian distribution $P_{lens}(r_{s},r_{200})$ to be combined with
the _MG-MAMPOSSt_ likelihood, following the approach of Ref. Pizzuti et al.
(2019). The additional distribution is assumed to be derived by an unbiased
(lensing) reconstruction, i.e. it is centered on the true values
$r_{200}^{t}=2.00\,\mbox{\rm Mpc},\,r_{s}^{t}=0.30\,\mbox{\rm Mpc}$ of the
cluster mass profile’s parameter. The standard deviations are defined assuming
the above average uncertainties from real lensing analyses:
$\sigma_{r_{200}}/r_{200}=0.1$ and $\sigma_{r_{s}}/r_{s}=0.3$. As for the
correlation, we assume $\rho=0.5$, checking that our results show only
negligible effects when varying $\rho$ over a reasonable range of values.

In BH theories , the lensing mass profile is determined by both parameters
$Y_{1}$ and $Y_{2}$. Indeed, as discussed in Section 2, the lensing mass
profile is given by a combination of $\Phi,\,\Psi$, Eq. (8). If $Y_{1}$ and
$Y_{2}$ are not equal to zero, a lensing mass reconstruction assuming standard
GR NFW model (which depends only on $r_{s}$ and $r_{200}$) leads to a biased
values of mass and concentration due to the modification of gravity. More in
detail, our simulated lensing distribution is written in terms of _effective
parameters_ :

$P_{\text{lens}}=P_{\text{lens}}(r_{200}^{e},r_{s}^{e}),$ (25)

where $r_{200}^{e},\,r_{s}^{e}$ are implicit functions of
$r_{200},r_{s},Y_{1},Y_{2}$. The relations between the two sets of parameters
are given by requiring that the values of the mass profiles
$M_{\text{NFW}}(r)$ and $M_{lens}(r)$, computed at each pair of radii
($r_{200}^{e},\,r_{s}^{e}$) should coincide,

$(r_{200}^{e})^{3}=\frac{r_{200}^{3}}{\log(1+c)-c/(1+c)}\left[\mathcal{G}_{1}(r_{200}^{e},r_{s}),\mathcal{G}_{2}(r_{200}^{e},r_{s},Y_{1},Y_{2})\right]\,,$
(26) $\frac{(r_{200}^{e})^{3}\left(\log
2-\frac{1}{2}\right)}{\log(1+c_{e})-c_{e}/(1+c_{e})}=.$
$\frac{r_{200}^{3}}{\log(1+c)-c/(1+c)}\left[\mathcal{G}_{1}(r_{s}^{e},r_{s}),\mathcal{G}_{2}(r_{s}^{e},r_{s},Y_{1},Y_{2})\right],$
(27)

where $c=r_{200}/r_{s}$ and $c_{e}=r_{200}^{e}/r_{s}^{e}$ and

$\displaystyle\mathcal{G}_{1}(r,r_{s})$ $\displaystyle=$
$\displaystyle\log\left(1+\frac{r}{r_{s}}\right)-\frac{r}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{-1}\,,$
(28) $\displaystyle\mathcal{G}_{2}(r,r_{s},Y_{1},Y_{2})$ $\displaystyle=$
$\displaystyle\frac{Y_{1}\left(\frac{r_{s}}{r}-1\right)-5Y_{2}\left(\frac{r_{s}}{r}+1\right)}{8(1+r_{s}/r)^{3}}.$
(29)

We carry out a Monte-Carlo-Markov-Chain sampling over $3\cdot 10^{5}$ points
in the space ($r_{200},r_{s},\beta_{\infty},$ $Y_{1},Y_{2}$) assuming flat
prior for each parameter. When $Y_{1},Y_{2}$ are considered, we excluded all
the combinations for which the resulting mass profiles become negative.
Moreover, $Y_{1}$ is assumed to be strictly larger than $-2/3$ to fulfill the
stability conditions highlighted in Saltas et al. (2018). The sampling has
been performed by using the joint distribution
$\log\mathcal{L}_{\text{tot}}=\log\mathcal{L}_{\text{dyn}}(r_{s},r_{200},\beta,Y_{1})+N\log
P_{\text{lens}}(r_{s}^{e},r_{200}^{e})$, where $N$ indicates the number of
clusters considered in the stacking analysis and
$\mathcal{L}_{\text{dyn}}(r_{s},r_{200},\beta,Y_{1})$ is the $N$ cluster-
combined _MG-MAMPOSSt_ likelihood. For each choice of
$(r_{200},r_{s},Y_{1},Y_{2})$, the lensing contribution to the total
likelihood is given by solving eqs. (27) to get $r_{200}^{e},\,r_{s}^{e}$ and
assigning the corresponding value of $P_{\text{lens}}(r_{s}^{e},r_{200}^{e})$.

The constraints at 95% C.L. on the BH parameters $Y_{1}$ and and $Y_{2}$ are
listed in Table 1 for different number of clusters considered in the combined
analysis. In fig. 7 we show the marginalized posterior distributions of
$Y_{1}$ and $Y_{2}$ in the case of a single halo. In this case, we obtain at
95% C.L.:

$Y_{1}\lesssim 4.18,\,\,\,\,\,\,\,\,\,\,Y_{2}=-0.22^{+0.69}_{-0.70}.$ (30)

Note that, as before, for $Y_{1}$ we consider only the upper limit, as the
lower bound is still given by the value settled by the stability conditions
$Y>-2/3$.

Figure 7:  Marginal distributions of $Y_{1}$ and $Y_{2}$ obtained from the lensing+dynamics likelihood for one cluster. In the lensing Gaussian we consider $\rho=0.5,\,\sigma_{r_{s}}=0.3\,r_{s}$ and $\sigma_{r_{200}}=0.1\,r_{200}$.  | Beyond Hordenski |  $f(R)$ gravity
---|---|---
| 600 tracers | 100 tracers | 600 tracers | 100 tracers
$N$ clusters | $Y_{1}$ | $Y_{2}$ | $Y_{1}$ | $Y_{2}$ | $|f_{R0}|$ | $|f_{R0}|$
1 | $\lesssim 4.18$ | $-0.22^{+0.69}_{-0.70}$ | $\lesssim 5.40$ | $-0.22^{+0.73}_{-1.05}$ | – | –
5 | $\lesssim 3.57$ | $-0.05^{+0.47}_{-0.49}$ | $\lesssim 4.71$ | $-0.01^{+0.40}_{-0.59}$ | $\lesssim 3.37\times 10^{-5}$ | $\lesssim 5.13\times 10^{-5}$
10 | $\lesssim 2.30$ | $0.01^{+0.27}_{-0.30}$ | $\lesssim 4.12$ | $-0.02^{+0.34}_{-0.44}$ | $\lesssim 1.12\times 10^{-5}$ | $\lesssim 3.24\times 10^{-5}$
15 | $0.49^{+1.28}_{-0.99}$ | $0.00^{+0.24}_{-0.25}$ | $\lesssim 3.96$ | $0.00^{+0.36}_{-0.39}$ | $\lesssim 9.51\times 10^{-6}$ | $\lesssim 2.43\times 10^{-5}$
20 | $0.52^{+1.07}_{-0.94}$ | $0.01^{+0.12}_{-0.15}$ | $\lesssim 3.73$ | $-0.01^{+0.30}_{-0.33}$ | $\lesssim 7.11\times 10^{-6}$ | $\lesssim 1.79\times 10^{-5}$
Table 1:  Constraints at 95% C.L forcasted from the _MG-MAMPOSSt_ method for
the parameters of the two modified gravity models presented in this paper. We
assumed an additional lensing distribution in the analysis modelled as a
Gaussian with $\rho=0.5,\,\sigma_{r_{s}}=0.3\,r_{s}$ and
$\sigma_{r_{200}}=0.1\,r_{200}$. From column two to column five: Beyond
Horndeski gravity. Column six and seven: chameleon $f(R)$ gravity. For both
models we show the results when using $\sim 600$ and $\sim 100$ cluster
members in the _MG-MAMPOSSt_ fit.

As expected, the results are consistent with the GR predictions and of the
same order of current cosmological bounds on those parameters derived with
similar methods. In particular, Ref. Sakstein et al. (2016) obtained
constraints of order unity by combining the stacked X-ray surface brightness
profiles from the XMM Cluster Surveyand weak lensing profiles from CFHTLenS of
58 clusters at redshift $0.1<z<1.2$. A much more stringent constraint has been
found by Salzano et al. (2017a) which however worked on a quite different
class of models where the potential $\Phi$ and $\Psi$ are modified in terms of
one single free parameter $\mathcal{Y}$ (although is still true that
$\Phi\neq\Psi$). By combining X-rays mass profiles and strong+weak lensing
analyses for 18 clusters studied within the CLASH collaboration, they found a
very tight bound of $\mathcal{Y}\lesssim 0.16$ at $2\,\sigma$. In view of the
above results, we remark the following:

  * •

Our forecasts are obtained by assuming an ideal situation with perfect control
of systematics for both lensing and dynamics analysis.

  * •

The uncertainties considered here are reliable, given the current available
data-sets, but slightly optimistic. Indeed, high-precision spectroscopic
measurements for $\sim 600$ cluster members could be obtained only for a small
fraction of observed galaxy clusters. Therefore, we consider a more realistic
situation of $\sim 100$ tracers in the _MAMPOSSt_ fit. The results are listed
in the fourth and fifth column of Table 1. As expected, in this case the
constraints are much weaker, especially for the parameter $Y_{1}$ even with a
large number of halos $N\gtrsim 15$.

  * •

The test carried out shows that galaxy cluster mass profiles are not able to
provide tight constraints for the Newtonian potential $\Phi$ in Beyond
Horndeski gravity (related to the coupling $Y_{1}$), compared to current
astrophysical constraints. Nevertheless, the combination of dynamics (galaxies
and gas)+lensing analyses in clusters promises to provide very useful
information on the parameter $Y_{2}$ in beyond Horndeski theories, for which
the only constraint obtained at cosmological level is the one derived in
Sakstein et al. (2016). Indeed, given the modification of the potentials Eq.
(6) and Eq. (7), the effective lensing mass for a NFW density profile reads
as,

$M_{\text{lens}}=M_{\text{NFW}}+\frac{r^{2}M_{200}\left[Y_{1}(r_{s}-r)-5Y_{2}(r_{s}+r)\right]}{4[\log(1+c)-c/(1+c)]}\frac{1}{(r_{s}+r)^{3}}$

$\equiv M_{\text{dyn}}+M_{2},$ (31)

with

$M_{2}=-\frac{5r^{2}M_{200}}{4[\log(1+c)-c/(1+c)]}\frac{Y_{2}(r_{s}+r)}{(r_{s}+r)^{3}}.$
(32)

We can see that the contribution of $Y_{2}$ to departures from GR is five time
stronger than $Y_{1}$. Moreover, the last term in Eq. (31) is always different
from zero for any finite $r$, as shown by the relative ratios of the lensing
mass profile in Fig. 8. Thus, we expect that joint Strong+Weak lensing and
dynamics analyses can constrain much better $Y_{2}$ than $Y_{1}$.

Figure 8:  Relative ratios for the effective lensing mass profiles
$(M_{lens}-M_{\text{NFW}})/M_{\text{NFW}}$ obtained in Beyond Horndeski
gravity with $Y_{1}=0.5,\,Y_{2}=0$ (blue) and $Y_{1}=0,\,Y_{2}=0.5$ (red). The
NFW profile parameters are $r_{200}=2.0\,\text{Mpc}$ and
$r_{s}=0.3\,\text{Mpc}$.

In this perspective, we further outline that:

  * •

The method we used to reconstruct the parameters from the simulated lensing
mass profile is surely not the same as performing a fit of the model to
convergence and shear maps. Our simplified procedure may underestimate the
uncertainties on the parameters.

  * •

In addition, the distribution $P_{\text{lens}}$ is estimated on the basis of
high precision joint Strong+Weak lensing analyses - like those performed by
Refs. Caminha et al. (2017); Umetsu et al. (2016) \- which can determine the
mass profile down to the cluster core. The structure of Eqs. (7),(6) and Eq.
(31) suggests that the knowledge of the inner shape of the profile provides
significant information in breaking the degeneracy between the model
parameters.

##  5 Application II: Chameleon screening

In this section we present the results obtained from the analysis of the
sample of synthetic halos when assuming a fifth force mediated by a chameleon
field $\phi$ in the _MG-MAMPOSSt_ method. In the general case, assuming that
the interior solution for the scalar field is close to zero up to the
screening radius, the magnitude of the additional force is totally determined
by the background field value $\phi_{\infty}$ and by the coupling to matter,
$\mathcal{Q}$. We first analyze our sample of synthetic phase spaces in this
framework, using Eq. (4) to model the chameleon profile. We perform a MCMC
sampling of the dynamics likelihood $\mathcal{L}_{\text{dyn}}$ over the full
parameter space $(r_{200},r_{s},\beta,\phi_{\infty},\mathcal{Q})$. As before,
$\beta\equiv(1-\beta_{\infty})^{-1/2}$, where $\beta_{\infty}$ is the
normalization of the Tiret anisotropy profile, Eq. (11).

In principle, while the chameleon screening mechanism recovers GR at small
scales, modifications of gravity are further Yukawa-suppressed for
$r\gtrsim\lambda_{\phi}\equiv
m^{-1}_{\phi}\sim\left[\mathcal{Q}R_{0}M_{p}/\phi_{\infty}\right]^{-1/2}$,
where $R_{0}$ is the background Ricci scalar. In order to satisfy local test
of gravity at Solar System scales, it has been shown that the background
interaction range $\lambda_{\phi}$ should be of the order of few Mpc. In this
simple treatment we neglect the Yukawa suppression of the chameleon field as
in the _MG-MAMPOSSt_ analysis we consider only galaxies lying in a projected
radius $R\sim r_{200}\sim\mathcal{O}(\text{Mpc})$. However, since in the
Jean’s equation the integral of the velocity dispersion extends far beyond the
viral radius, for large $\mathcal{Q}$ the motion of member galaxies could be
affected by the contribution of the tails of the integral even if the cluster
is totally screened (i.e. $S\simeq r_{200}$). In order to avoid this problem,
we make a conservative assumption that the effect of the chameleon field is
zero everywhere if $S>1.5\,r_{200}$. A complete solution of the field equation
will be investigated elsewhere.

Since in chameleon gravity the lensing mass profile is unaffected by the
additional fifth force, the contribution of the lensing probability
distribution $P_{L}(r_{s},\,r_{200})$ affects only the NFW mass profile
parameters $r_{200}$, $r_{s}$. We thus re-sampled our parameter space over the
joint (log) Likelihood
$\log[\mathcal{L}_{\text{dyn}}(r_{s},r_{200},\beta,\phi_{\infty},\mathcal{Q})]+\log[P_{L}(r_{s},\,r_{200})]$.
For our reference analysis, as before, we consider a massive halo with
$r_{200}=2.0\,\text{Mpc},\,r_{s}=0.3\,\text{Mpc}$ and $\sim 600$ tracers in
the _MG-MAMPOSSt_ fit and
$\sigma_{r_{200}}/r_{200}=0.1,\,\sigma_{r_{s}}/r_{s}=0.3,\,\rho=0.5$ in the
additional (lensing) distribution. The result of the MCMC for the MG
parameters are shown in the upper plots of Fig. 9.

Figure 9:  Top: results for the MCMC run in chameleon gravity for one halo in
the sample. The shaded areas show the allowed regions in the plane
$(\mathcal{Q}_{2}=\mathcal{Q}/(1+\mathcal{Q}),\phi_{2}=1-\exp(\phi_{\infty}M_{p}/10^{-4}))$
at $3\sigma$ (outer region) and $2\sigma$ (innermost region). Left panel:
results from the _MG-MAMPOSSt_ Likelihood. Right panel: results from the joint
lensing+dynamics analysis with $\sigma r_{200}/r_{200}=0.1$, $\sigma
r_{s}/r_{s}=0.3,\,\rho=0.5$ in the additional lensing distribution. Bottom:
joint lensing+dynamics over the combination of 10 halos. In this and the
following figures we use red-yellow for _MG-MAMPOSSt_ alone, blue-light blue
for the joint lensing+dynamics.

In order to explore the full parameter range, we plot the rescaled variables
$\mathcal{Q}_{2}=\mathcal{Q}/(1+\mathcal{Q}),\,\,\phi_{2}=1-\exp(\phi_{\infty}M_{p}/10^{-4})$,
which run in the interval $[0,1)$. The results are similar to e.g. Refs.
Terukina et al. (2014b); Wilcox et al. (2015b), which performed joint lensing
and X-rays analyses of 58 clusters and of the massive Coma cluster
respectively. For small coupling factors, all field values are allowed by
galaxy cluster dynamics with the screening mechanism active at small
$\phi_{\infty}$. As the screening radius depends also on $\mathcal{Q}$, at
large $\mathcal{Q}$ the chameleon screening becomes very efficient also for
large field values. In the intermediate region, however, $\phi_{2}\to 1$ is no
more allowed if $\mathcal{Q}\gtrsim 0.6$ as the cluster will be only partially
screened and the modification of gravity cannot be entirely suppressed.
Already with a single cluster, the _MG-MAMPOSSt_ method alone and the
combination of _MG-MAMPOSSt_ with additional lensing distribution provide
useful information to investigate this intermediate range of values.

It is worth to point out that the excluded region found with _MG-MAMPOSSt_ is
shifted towards larger $\mathcal{Q}_{2}$ w.r.t. what found by Refs. Terukina
et al. (2014b); Wilcox et al. (2015b). This is not surprising as X-rays
analyses and Jeans’ equation analyses perceive the effect of gravity in
different ways, despite they are sensitive to the same gravitational
potential. In particular, as mentioned before, the radial velocity dispersion
$\sigma^{2}_{r}(r)$ is affected by the mass profile in the outskirt of the
cluster (i.e. the integral of Eq. (12) used in the _MG-MAMPOSSt_ method
extends from $r$ to $\infty$). On the other side, the solution of hydrostatic
equilibrium for the electron gas pressure at a given radius is given by (e.g.
Ref. Wilcox et al. (2015b)):

$P_{e}(r)=P_{e,0}+\mu
m_{p}\int_{0}^{r}{n_{e}(r)\left[-\frac{GM(<r)}{r^{2}}-\frac{\mathcal{Q}}{M_{p}}\frac{d\phi(r)}{dr}\right]},$
(33)

where $n_{e}(r)$ is the electron density profile, $m_{p}$ the proton mass
$P_{e,0}$ the gas pressure at $r=0$ and $\mu$ the mean molecular weight. Thus,
in this case the integral is from $0$ to $r$, i.e. the pressure is determined
by the contribution _up to_ a given radius.

Our results illustrate that the _MG-MAMPOSSt_ procedure, together with
strong+weak lensing analyses, is an interesting method to test chameleon
gravity at cluster scales – with present and future observations, the
combination of different sets of data and methodologies (i.e. X-rays mass
profile determinations) can improve the constraints on the allowed region in
the plane $(\mathcal{Q},\phi_{\infty})$ and thus on the viability of those
theories at cosmological level. As for Beyond Horndeski model, we study how
our results change when combining the information of several clusters, by
performing a MCMC sampling of
$\log\mathcal{L}_{tot}=\log\mathcal{L}_{\text{dyn}}(r_{s},r_{200},\beta,\phi_{\infty},\mathcal{Q})+N\log
P_{lens}(r_{s},r_{200})$, where $N$ is the number of halos considered in the
stack. As shown in the bottom plot of Fig. 9, the analysis of 10 clusters
produce an improvement of the constraints for $\mathcal{Q}_{2}\lesssim 0.4$,
tightening the allowed region in the parameter space.

Figure 10:  $3\sigma$ (outer region) and $2\sigma$ (innermost region) regions
for the analysis of one cluster with $N=100$ (top plots) and $N=50$ (bottom
plots) tracers considered in the _MG-MAMPOSSt_ fit.The color code is the same
as in Fig. 9, i.e. red-yellow for _MG-MAMPOSSt_ alone, blue-light blue for the
joint lensing+dynamics.

We point out again that the availability of 600 galaxies to be used in the
_MG-MAMPOSSt_ will be limited to few clusters from present and future surveys.
We thus reduce the number of tracers $N$ from 600 to 100 and we investigate
the behavior of the iso-probability contours in this case, shown in the upper
plots of Fig. 10. While the _MG-MAMPOSSt_ analysis alone produces weaker
constraints with respect to the reference case $N=600$, as expected when
decreasing the statistics, the joint lensing+dynamics contours are almost
identical to the previous analysis. This indicates that the constraining power
is mostly related to the additional information provided on the mass profile
parameters $r_{200},r_{s}$, confirming what already found in Appendix A of
Ref. Pizzuti et al. (2020) in the simple case of linear $f(R)$. Keeping fixed
$P_{\text{lens}}(r_{s},\,r_{200})$, we generate a new phase space considering
only 50 tracers to push further this ideal exercise. As we can see in the left
bottom plot of Fig. 10, the statistics is insufficient to provide bounds in
$(\mathcal{Q},\phi_{\infty})$ with the MG-MAMPOSSt method. Now, the weakening
of the constraints affects also the joint lensing+dynamics distribution (right
plot). Although useful information can still be extracted in this limiting
situation (e.g. by considering several clusters and suitably combining the
member galaxy dynamics with real lensing and X-rays data), it is worth to
notice that we are neglecting all the systematic effects. Indeed, such a small
number of galaxies for each cluster cannot be enough to correctly track the
underlying total gravitational potential – incompleteness of the member
galaxies sample as well as the effect of interlopers (i.e. galaxies in the
field of view which do not belong to the cluster) can produce a bias in
dynamic mass estimation if not properly taken into account (see Ref. Mamon et
al. (2013)).

Whereas for the beyond Horndeski model implemented in _MG-MAMPOSSt_ the new
degrees of freedom $Y_{1}\,Y_{2}$ are scale-independent, this is no more true
for chameleon gravity. The screening mechanism depends on the structure of the
halo density profile, as shown by Eq. (5) – thus, the analysis of the member
galaxy dynamics for clusters with different $r_{s},\,\,r_{200}$ should produce
different constraints. For this reason, we consider four halos generated with
other combination of $r_{s},\,r_{200}$. The new phase spaces are in two pairs
characterized by a concentrations $c=8.7$ and $c=5$ respectively, but obtained
by changing $r_{200}$ for one halo in the pair and $r_{s}$ for the other halo
with respect to the reference case
$r_{200}=2.0\,\text{Mpc},\,\,r_{s}=0.3\,\text{Mpc},\,\,c=6.67$. We apply again
the _MG-MAMPOSSt_ technique on phase spaces with $\sim 600$ particles within
$r_{200}$ drawn from those halos – we derive the marginalized iso-probability
contours for the modified gravity parameters, shown in Fig. 11. As expected,
changing both $r_{s}$, $r_{200}$ has an impact on the constrained allowed
region, with the strongest effect found for the smallest value of the
concentration and the smallest virial radius. In this case the _MG-MAMPOSSt_
analysis alone provides almost the same constraints as the joint
lensing+dynamics distribution obtained for the reference case (right plot of
Fig. 9).

Figure 11:  Result of the _MG-MAMPOSSt_ analysis (dynamics only) for single
halos with different values of $r_{s}$ and $r_{200}$. Color-codes is the same
as previous plots. The left upper and left bottom plots refers to a
concentration $c=5$ obtained changing $r_{200}$ and $r_{s}$ respectively from
the reference case. The upper right and bottom right plots show the result for
a higher concentration parameter $c=8.7$. Given the same concentration, the
stringent bounds are obtained for the less massive halo, i.e. the smaller
$r_{200}$.

####  $f(R)$ gravity

We now restrict to $f(R)$ gravity, a sub-class of chameleon theories where the
coupling parameter is fixed to $\mathcal{Q}=1/\sqrt{6}$. Constraining this
model translates into bounds on the background value of the scalaron field
$|f_{R0}|$, which mediates the additional fifth force and it is related to
$\phi_{\infty}$ as

$\phi_{\infty}=-\sqrt{\frac{3}{2}}\log(1+f_{R0}),$ (34)

in units of $M_{p}c^{2}$.

For our reference analysis, we consider again a phase space of $\sim 600$
particles within $r_{200}$ from an halo with $r_{200}=2.0\,\text{Mpc}$,
$r_{s}=0.3\,\text{Mpc}$. In the simulated lensing Gaussian distribution we
choose the same values of $\sigma_{r_{200}}$, $\sigma_{r_{s}}$ and $\rho$ as
before. We then vary the number of halos considered in the MCMC run, the
lensing standard deviations of $r_{200}$ and $r_{s}$ and the number of tracers
in the dynamics analysis.

When considering a single halo, we found that $\phi_{\infty}$ is unconstrained
at $95\%$ C.L. in both cases of $\sim 600$ and $\sim 100$ tracers in projected
phase space and for $\sigma_{r_{200}}/r_{200}=0.1,\,0.07$. Moreover, changing
$\sigma_{r_{s}}$ in the lensing distribution produces a negligible effect on
the joint likelihood. This is due to the degeneracy between $r_{200}$ and the
chameleon field,confirming what found in the simplified picture of linear
$f(R)$ by Ref. Pizzuti et al. (2020). Indeed, if we fix the coupling constant
$\mathcal{Q}=1/\sqrt{(}6)$, very large $\phi_{\infty}$ produces a maximum
enhancement of the GR gravitational potential by a factor $1/3$ which can be
compensated by slightly decreasing the overall radius of the halo. More in
particular, a cluster with a smaller $r_{200}$ in a strong modified $f(R)$
gravity ($\phi_{\infty}\gtrsim 10^{-4}$) produce the same effect of a cluster
with a larger $r_{200}$ in a GR scenario. In order to obtain a (weak) bound on
$\phi_{\infty}$, the uncertainties on $r_{200}$ in the lensing Gaussian should
be reduced down to $5\%$. In this case we found $\phi_{\infty}\lesssim
9.52\times 10^{-5}$ at $\Delta\chi^{2}=4.0$, corresponding to
$|f_{R0}|\lesssim 7.77\times 10^{-5}$.

As for the general chameleon model, we perform a stacked analysis to study the
constraining power of our method in the ideal situation. In Fig. 12 we plot
the marginalized likelihood obtained by Monte-Carlo sampling over a combined
lensing+dynamics likelihood of 10 halos. The black-dashed, blue and red lines
indicate different values of $\sigma_{r_{200}}/r_{200}$ in the lensing
Gaussian, 10% (black line) and 5% (red line). The uncertainties on $r_{s}$ are
fixed to 30% and the correlation to $\rho=0.5$. The brown curve is obtained
with 100 tracers considered in the _MG-MAMPOSSt_ fit in each of the 10 halos.

For $\sigma_{r_{200}}=0.1\,r_{200}$ and $\sigma_{r_{200}}=0.05\,r_{200}$ the
combined distributions are almost identical, indicating that the information
which can be obtained with our method for this particular model is already
saturated with $\sim 10$ galaxy clusters. This extends the results in Appendix
A of Ref. Pizzuti et al. (2020) to the case of chameleon $f(R)$ gravity. For
the reference analysis ($\sigma_{r_{200}}/r_{200}=0.1$ in the lensing
distribution, red dotted line in of Fig. 12) and 10 halos considered in the
combined likelihood we obtain (in terms of $f_{R0}$):

$\displaystyle|f_{R0}|$ $\displaystyle\leq$ $\displaystyle 8.26\times
10^{-6}\,\,\,\,1\sigma,$ (35) $\displaystyle|f_{R0}|$ $\displaystyle\leq$
$\displaystyle 1.12\times 10^{-5}\,\,\,\,2\sigma.$ (36)

When decreasing the number of tracers down to 100, the constraints become
slightly weaker, $|f_{R0}|\leq 3.24\times 10^{-5}$ at $2\sigma$, but still
good if compared with the current cosmological constraints on $f(R)$. If the
number of halos in the combined likelihood is extended up to 20, the
constraints improve to $|f_{R0}|\leq 1.79\times 10^{-5}$ and $|f_{R0}|\leq
7.11\times 10^{-6}$ at $2\sigma$ for 100 tracers and 600 tracers in the _MG-
MAMPOSSt_ fit, respectively. The results obtained varying the number of
clusters considered in the analysis are listed in the last two columns of
Table 1.

As mentioned before, the bounds can be substantially tightened by considering
a less efficient screening mechanism, i.e. a smaller-size halo. The blue curve
in Fig. 12 is obtained for 10 halos with $r_{200}=1.5\,\mbox{\rm Mpc}$ and
$r_{s}=0.3\,\mbox{\rm Mpc}$. In this case, the upper limits become

$\displaystyle|f_{R0}|$ $\displaystyle\leq$ $\displaystyle 3.48\times
10^{-6}\,\,\,\,1\sigma,$ (37) $\displaystyle|f_{R0}|$ $\displaystyle\leq$
$\displaystyle 5.40\times 10^{-6}\,\,\,\,2\sigma.$ (38)

The constraints further strengthen to $|f_{R0}|\leq 3.56\times 10^{-6}$ at
$2\sigma$ in the case of 20 halos. The bounds achieved by this simple exercise
with a reasonable number of clusters are of the same order of magnitude of the
results obtained from other cosmological probes (see e.g. Wilcox et al.
(2015b); Cataneo et al. (2016a); Jana & Mohanty (2019)) for the specific class
of $f(R)$ models. As already mentioned, it will be interesting to further
evaluate the effect of X-rays analyses which could provide additional
information to the dynamical parameters. We point out again that our results
are obtained neglecting all systematics, as a proof of concept of the
statistical power of the _MG-MAMPOSSt_ method.

Figure 12:  Marginalized distribution of the chameleon field in the case of
$f(R)$ gravity ($\mathcal{Q}=1/\sqrt{6}$) obtained from the combined, joint
lensing+dynamics likelihood for 10 clusters. In the lensing Gaussian we
consider $\rho=0.5,\,\sigma_{r_{s}}=0.3\,r_{s}$ and a variable
$\sigma_{r_{200}}$. Black dotted line refers to the reference analysis with
600 tracers in _MG-MAMPOSSt_ and $\sigma_{r_{200}}/r_{200}=0.1$. The brown
curve shows the case of 100 tracers while the blue distribution is for a phase
space of 600 tracers but derived from a cluster with $r_{200}=1.5\,\mbox{\rm
Mpc}$. Finally, the red line refers to $\sigma_{r_{200}}/r_{200}=0.05$.

##  6 Summary and Conclusions

We have presented a new extension of the _MG-MAMPOSSt_ code, aimed at
constraining modified gravity models with the dynamics of the member galaxies
in clusters. In particular, we have implemented two parametrisations for the
Newtonian potential $\Phi$, which governs the dynamics of member galaxies,
based on Beyond Hordenski and generic chameleon scalar-tensor theories and NFW
density profile for matter perturbations. These models cover the range of
viable scalar-tensor theories after the recent measurement of the speed of
gravitational waves. The code is based on the original _MAMPOSSt_ method
developed by Ref. Mamon et al. (2013) which applies a Maximum Likelihood
approach by solving the spherical Jeans’ equation in the projected phase space
of member galaxies $(R,v_{z})$. As shown in Section 3, _MG-MAMPOSSt_ can be
used to fit together the mass and number density profile parameters, the
velocity anisotropy parameter plus the additional degrees of freedom related
to the modified gravity theory under study.

The code has been equipped with a simple MCMC module based on a Metropolis-
Hastings algorithm for an efficient exploration of the parameter space of the
theories. A run takes few hours to compute the likelihood for $\sim 10^{5}$
sampled values of the parameters for the case of all the implemented models.

Under the assumptions of dynamical relaxation and spherical symmetry, we have
illustrated the code’s capabilities for a set of synthetic isolated Dark
Matter halos, generated by using the ClusterGEN code of Ref. Pizzuti et al.
(2019). Working within this ideal setup, in our analysis we investigated the
statistical degeneracy between the model parameters for Beyond Horndeski and
Chameleon gravity that could obscure future realistic constraints in this
context. We have further produced forecasts on those parameters in view of
current and future observational surveys. The results are summarized as
follows:

In the case of Beyond Horndeski gravity we found that the resulting likelihood
computed with _MG-MAMPOSSt_ based on the dynamics of member galaxies, exhibits
a considerable degeneracy between the NFW parameters and the coupling constant
$Y_{1}$. An important upshot in this regard is that additional information is
required in order to break this degeneracy and constrain this class of
modified gravity models.

Since for this family of models lensing is also affected as compared with GR,
we simulated additional information on the cluster’s effective mass profile,
using reliable strong+weak lensing information from the literature. We showed
that the forecasted constraints obtained on the fifth force coupling $Y_{1}$
(associated to the Newtonian potential $\Phi$) are order unity, significantly
weaker to other astrophysical probes. In particular, by considering 20 halos’
phase spaces with 600 tracers in the _MG-MAMPOSSt_ fit plus a lensing Gaussian
distribution with $\rho=0.5,\,\sigma_{r_{s}}=0.3\,r_{s}$ and
$\sigma_{r_{200}}=0.1\,r_{200}$, we obtain $Y_{1}=0.52^{+1.07}_{-0.94}$. When
reducing the number of galaxies in each phase space to 100 the constraints
weaken to $Y_{1}\lesssim 3.73$. However, as we further showed, the joint
dynamics+lensing analysis can be used to place tight constraints on the
parameter $Y_{2}$ which appears in the relativistic potential $\Psi$ and it
has a stronger impact to the observed mass profile with respect to $Y_{1}$. In
this case the aforementioned analysis we found $Y_{2}=-0.01^{+0.12}_{-0.15}$
and $Y_{2}=-0.01^{+0.30}_{-0.33}$ for 600 and 100 tracers in the fit
respectively.

For chameleon theories we considered the general case with two free
parameters, the fifth-force coupling strength $\mathcal{Q}$ and the background
value of the scalar field $\phi_{\infty}$. The _MG-MAMPOSSt_ code, alone or
together with additional (lensing) information, provides useful insights on
the allowed region of the parameter space,
$\mathcal{Q}_{2}=\mathcal{Q}/(1+\mathcal{Q}),\phi_{2}=1-\exp(\phi_{\infty}M_{p}/10^{-4})$.
We found that, given the dependence of the velocity dispersion on the
gravitational potential, the module implemented in _MG-MAMPOSSt_ can be used
to get complementary information to X-rays mass reconstructions, even if they
are both sensitive to the time-time potential $\Phi$.

In this class of models, lensing determinations are not affected by the fifth
force, thus they can be used as a prior information on the cluster parameters
$r_{s}\,,r_{200}$. Those priors have a strong impact on the constraining power
of the method. Indeed, we found that decreasing the number of tracers from 600
to 100 in the _MG-MAMPOSSt_ fit but keeping the same simulated lensing
distribution produces a negligible effect on the two dimensional iso-
probability contours in the plane ($\mathcal{Q}_{2}$, $\phi_{2}$).

For our reference analysis of a phase space with 600 galaxies and
$\rho=0.5,\,\sigma_{r_{s}}=0.3\,r_{s}$, $\sigma_{r_{200}}=0.1\,r_{200}$ in the
lensing Gaussian we found that the region corresponding to
$0.45\lesssim\mathcal{Q}_{2}\lesssim 0.7$ and $\phi_{2}\gtrsim 0.6$ is
excluded at more than $3\sigma$. The combined likelihood of 10 halos increases
the range of the prohibited area up to $\phi_{2}\gtrsim 0.3$ for
$\mathcal{Q}_{2}$ between $0.35$ and $0.6$.

The chameleon screening depends on the matter density perturbations, producing
a stronger damping of the fifth force in more massive clusters. We considered
four different NFW halos characterized by distinct values of $r_{s}$ and
$r_{200}$. We confirmed that smaller-size, less concentrated halos provide
better constraints due to a less-efficient screening mechanism. More in
detail, for a cluster with $r_{s}=0.3\,\mbox{\rm Mpc}$ and
$r_{200}=1.5\,\mbox{\rm Mpc}$ the _MG-MAMPOSSt_ fit alone provide almost the
same constraints as the joint lensing+dynamics analysis of a cluster with
$r_{s}=0.3\,\mbox{\rm Mpc}$ and $r_{200}=2.0\,\mbox{\rm Mpc}$.

As a final step, we have focused on the specific sub-class of chameleon $f(R)$
gravity (i.e. $\mathcal{Q}=1/\sqrt{6}$,
$\phi_{\infty}=-\sqrt{\frac{3}{2}}\log(1+f_{R0})$), in order to forecast the
value of the background scalaron field. We obtain, for the reference case,
$|f_{R0}|\lesssim 1.12\times 10^{-5}$ and $|f_{R0}|\lesssim 7.11\times
10^{-6}$ at 95% C.L. from the joint lensing+dynamics analysis of 10 clusters
and 20 clusters, respectively. The bounds can be improved when using less
massive clusters – in particular we found $|f_{R0}|\leq 5.40\times 10^{-6}$
and $|f_{R0}|\leq 3.56\times 10^{-6}$ for 10 and 20 halos with
$r_{s}=0.3\,\mbox{\rm Mpc}$ and $r_{200}=1.5\,\mbox{\rm Mpc}$.

In conclusion, _MG-MAMPOSSt_ offers a tool to probe GR and viable scalar-
tensor theories at the scale of a galaxy clusters. A key aspect that has to be
accounted for is the statistical degeneracy between the mass profile
parameters and the MG parameters, at least for the models currently
implemented in the code. As we discussed earlier, the breaking of degeneracy
requires combination of additional probes such as lensing or X-ray analyses,
in order to provide tight and accurate constraints. Both ground-based (e.g.
LSST) and space telescopes (e.g. Euclid), together with next generation
spectrographs will provide a large amount of high quality imaging and
spectroscopic data for several galaxy clusters, allowing for a wide range of
applications of our method.

The code is still in development: we aim at extending the parametrisation of
the gravitational potential to other mass density profiles, which may better
describe the total matter distribution in galaxy clusters, as well as to other
viable modified gravity models. It is worth pointing out that the results
presented here are obtained assuming the ideal conditions required by the
limitations of our method. When dealing with real data, several systematics
can affect the measurement producing spurious detection of modified gravity.
For example, as shown in Ref. Pizzuti et al. (2020), lack of dynamical
relaxation and deviation from spherical symmetry have a relevant impact on the
projected phase space and should be taken into account by adopting suitable
selection criteria. The application of the _MG-MAMPOSSt_ method on the
observed phase space of relaxed galaxy clusters will be performed in a
forthcoming paper.

## Acknowledgements

LP is partially supported by a 2019 "Research and Education" grant from
Fondazione CRT. The OAVdA is managed by the Fondazione Cleŕment Fillietroz-
ONLUS, which is supported by the Regional Government of the Aosta Valley, the
Town Municipality of Nus and the "United́es Communes valdotaines Mont-Eḿilius.
I.D.S. is funded by European Structural and Investment Funds and the Czech
Ministry of Education, Youth and Sports (Project CoGraDS —
CZ.02.1.01/0.0/0.0/15_003/0000437).

## References

  * Amendola et al. (2013) Amendola L., Kunz M., Motta M., Saltas I. D., Sawicki I., 2013, Phys. Rev. D, 87, 023501
  * Amendola et al. (2018) Amendola L., Kunz M., Saltas I. D., Sawicki I., 2018, Phys. Rev. Lett., 120, 131101
  * Amendola et al. (2020) Amendola L., Bettoni D., Pinho A. M., Casas S., 2020, Universe, 6, 20
  * Babichev & Lehébel (2018) Babichev E., Lehébel A., 2018, JCAP, 1812, 027
  * Babichev et al. (2016) Babichev E., Koyama K., Langlois D., Saito R., Sakstein J., 2016, Class. Quant. Grav., 33, 235014
  * Baker et al. (2017) Baker T., Bellini E., Ferreira P. G., Lagos M., Noller J., Sawicki I., 2017
  * Ben Achour et al. (2016) Ben Achour J., Crisostomi M., Koyama K., Langlois D., Noui K., Tasinato G., 2016, JHEP, 12, 100
  * Biviano & Salucci (2006) Biviano A., Salucci P., 2006, A&A, 452, 75
  * Biviano et al. (2013) Biviano A., et al., 2013, A&A, 558, A1
  * Budzynski et al. (2012) Budzynski J. M., Koposov S. E., McCarthy I. G., McGee S. L., Belokurov V., 2012, MNRAS, 423, 104
  * Burrage & Sakstein (2016) Burrage C., Sakstein J., 2016, JCAP, 11, 045
  * Caminha et al. (2017) Caminha G. B., et al., 2017, A&A, 607, A93
  * Cataneo et al. (2015) Cataneo M., et al., 2015, Physical Review D - Particles, Fields, Gravitation and Cosmology, 92
  * Cataneo et al. (2016a) Cataneo M., Rapetti D., Lombriser L., Li B., 2016a, preprint,  (arXiv:1607.08788)
  * Cataneo et al. (2016b) Cataneo M., Rapetti D., Lombriser L., Li B., 2016b, JCAP, 1612, 024
  * Creminelli & Vernizzi (2017) Creminelli P., Vernizzi F., 2017
  * Creminelli et al. (2018) Creminelli P., Lewandowski M., Tambalo G., Vernizzi F., 2018, JCAP, 1812, 025
  * Crisostomi & Koyama (2018) Crisostomi M., Koyama K., 2018, Phys. Rev., D97, 021301
  * Crisostomi et al. (2019) Crisostomi M., Lewandowski M., Vernizzi F., 2019
  * Dima & Vernizzi (2018) Dima A., Vernizzi F., 2018, Phys. Rev., D97, 101302
  * Ezquiaga & Zumalacárregui (2017) Ezquiaga J. M., Zumalacárregui M., 2017
  * Ferraro et al. (2011) Ferraro S., Schmidt F., Hu W., 2011, Phys. Rev. D, 83, 063503
  * Host et al. (2009) Host O., Hansen S. H., Piffaretti R., Morandi A., Ettori S., Kay S. T., Valdarnini R., 2009, ApJ, 690, 358
  * Ishak (2019) Ishak M., 2019, Living Rev. Rel., 22, 1
  * Jain et al. (2016) Jain R. K., Kouvaris C., Nielsen N. G., 2016, Phys. Rev. Lett., 116, 151103
  * Jana & Mohanty (2019) Jana S., Mohanty S., 2019, Physical Review D, 99
  * Kobayashi (2019) Kobayashi T., 2019
  * Kobayashi & Hiramatsu (2018) Kobayashi T., Hiramatsu T., 2018, Phys. Rev., D97, 104012
  * Kobayashi et al. (2015) Kobayashi T., Watanabe Y., Yamauchi D., 2015, Phys. Rev., D91, 064013
  * Koyama & Sakstein (2015) Koyama K., Sakstein J., 2015, Phys. Rev. D, 91, 124066
  * Langlois (2018) Langlois D., 2018
  * Langlois et al. (2017) Langlois D., Mancarella M., Noui K., Vernizzi F., 2017, JCAP, 1705, 033
  * Lombriser et al. (2012) Lombriser L., Slosar A., Seljak U., Hu W., 2012, Phys.Rev., D85, 124038
  * Mamon & Łokas (2005) Mamon G. A., Łokas E. L., 2005, MNRAS, 363, 705
  * Mamon et al. (2010) Mamon G. A., Biviano A., Murante G., 2010, A&A, 520, A30
  * Mamon et al. (2013) Mamon G. A., Biviano A., Boué G., 2013, MNRAS, 429, 3079
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Pizzuti et al. (2016) Pizzuti L., et al., 2016, JCAP, 1604, 023
  * Pizzuti et al. (2017) Pizzuti L., et al., 2017, JCAP, 1707, 023
  * Pizzuti et al. (2019) Pizzuti L., Saltas I. D., Casas S., Amendola L., Biviano A., 2019, Mon. Not. Roy. Astron. Soc., 486, 596
  * Pizzuti et al. (2020) Pizzuti L., Sartoris B., Borgani S., Biviano A., 2020, JCAP, 04, 024
  * Rapetti et al. (2010) Rapetti D., Allen S. W., Mantz A., Ebeling H., 2010, Monthly Notices of the Royal Astronomical Society, 406, 1796
  * Rapetti et al. (2011) Rapetti D., Allen S. W., Mantz A., Ebeling H., 2011, Progress of Theoretical Physics Supplement, 190, 179
  * Sakstein (2013) Sakstein J., 2013, Phys. Rev. D, 88, 124013
  * Sakstein (2015a) Sakstein J., 2015a, Phys. Rev. Lett., 115, 201101
  * Sakstein (2015b) Sakstein J., 2015b, Phys. Rev., D92, 124045
  * Sakstein (2018) Sakstein J., 2018, Int. J. Mod. Phys., D27, 1848008
  * Sakstein & Jain (2017) Sakstein J., Jain B., 2017
  * Sakstein et al. (2016) Sakstein J., Wilcox H., Bacon D., Koyama K., Nichol R. C., 2016, JCAP, 1607, 019
  * Saltas & Lopes (2019) Saltas I. D., Lopes I., 2019, Phys. Rev. Lett., 123, 091103
  * Saltas et al. (2018) Saltas I. D., Sawicki I., Lopes I., 2018, JCAP, 1805, 028
  * Salzano et al. (2016) Salzano V., Mota D. F., Dabrowski M. P., Capozziello S., 2016, JCAP, 1610, 033
  * Salzano et al. (2017a) Salzano V., Mota D. F., Capozziello S., Donahue M., 2017a, Physical Review D, 95
  * Salzano et al. (2017b) Salzano V., Mota D. F., Capozziello S., Donahue M., 2017b, Phys. Rev., D95, 044038
  * Schmidt et al. (2009) Schmidt F., Vikhlinin A., Hu W., 2009, Phys. Rev. D, 80, 083505
  * Stark et al. (2019) Stark A., Miller C. J., Halenka V., 2019, The Astrophysical Journal, 874, 33
  * Terukina & Yamamoto (2012a) Terukina A., Yamamoto K., 2012a, Phys. Rev. D, 86, 103503
  * Terukina & Yamamoto (2012b) Terukina A., Yamamoto K., 2012b, Phys. Rev. D, 86, 103503
  * Terukina et al. (2014a) Terukina A., Lombriser L., Yamamoto K., Bacon D., Koyama K., Nichol R. C., 2014a, JCAP, 1404, 013
  * Terukina et al. (2014b) Terukina A., Lombriser L., Yamamoto K., Bacon D., Koyama K., Nichol R. C., 2014b, Journal of Cosmology and Astroparticle Physics, 2014, 013–013
  * Tiret et al. (2007) Tiret O., Combes F., Angus G. W., Famaey B., Zhao H. S., 2007, A&A, 476, L1
  * Umetsu et al. (2016) Umetsu K., Zitrin A., Gruen D., Merten J., Donahue M., Postman M., 2016, ApJ, 821, 116
  * Wilcox et al. (2015a) Wilcox H., et al., 2015a, Mon. Not. Roy. Astron. Soc., 452, 1171
  * Wilcox et al. (2015b) Wilcox H., et al., 2015b, Monthly Notices of the Royal Astronomical Society, 452, 1171–1183
  * Wojtak et al. (2009) Wojtak R., Łokas E. L., Mamon G. A., Gottlöber S., 2009, Monthly Notices of the Royal Astronomical Society, 399, 812
  * Zumalacárregui & García-Bellido (2013) Zumalacárregui M., García-Bellido J., 2013