#  Testing generalized scalar-tensor theories of gravity with clusters of
galaxies

Balakrishna S. Haridasu<EMAIL_ADDRESS>SISSA-International School
for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy  INFN, Sezione di
Trieste, Via Valerio 2, I-34127 Trieste, Italy  IFPU, Institute for
Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy
Purnendu Karmakar<EMAIL_ADDRESS>Dipartimento di Fisica,
Sapienza Universitá di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  Marco De
Petris<EMAIL_ADDRESS>Dipartimento di Fisica, Sapienza
Universitá di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  INAF - Osservatorio
Astronomico di Roma, Via Frascati 33, 00040, Monteporzio Catone, Roma, Italy
INFN- Sezione di Roma 1, P.le Aldo Moro 2, 00185, Roma, Italy  Vincenzo F.
Cardone<EMAIL_ADDRESS>INAF - Osservatorio Astronomico di Roma, Via
Frascati 33, 00040, Monteporzio Catone, Roma, Italy  INFN- Sezione di Roma 1,
P.le Aldo Moro 2, 00185, Roma, Italy  Roberto Maoli
<EMAIL_ADDRESS>Dipartimento di Fisica, Sapienza Universitá di
Roma, P.le Aldo Moro 2, 00185, Roma, Italy  INAF - Osservatorio Astronomico di
Roma, Via Frascati 33, 00040, Monteporzio Catone, Roma, Italy

( at of Roma)

###### Abstract

We test the generalized scalar-tensor theory in static systems, namely galaxy
clusters. The Degenerate higher-order scalar-tensor (DHOST) theory modifies
the Newtonian potential through effective Newtonian constant and $\Xi_{1}$
parameter in the small scale, which modifies the hydrostatic equilibrium. We
utilize the well-compiled X-COP catalog consisting of 12 clusters with Intra
Cluster Medium (ICM) pressure profile by Sunyaev-Zel’dovich effect data and
temperature profile by X-ray data for each cluster. We perform a fully
Bayesian analysis modeling Navarro–Frenk–White (NFW) for the mass profile and
the simplified Vikhlinin model for the electron density. Carefully selecting
suitable clusters to present our results, we find a mild to moderate, i.e,
$\sim 2\sigma$ significance for a deviation from the standard scenario in 4 of
the clusters. However, in terms of Bayesian evidence, we find either
equivalent or mild preference for GR. We estimate a joint constraint of
$\Xi_{1}=-0.030\pm 0.043$ using 8 clusters, for a modification from a
$\Lambda$CDM scenario. This limit is in very good agreement with theoretical
ones and an order of magnitude more stringent than the previous constraint
obtained using clusters. We also quote a more conservative limit of
$\Xi_{1}=-0.061\pm 0.074$. Finally, we comment on the tentative redshift
dependence ($\Xi_{1}(z)$), finding a mild preference ($\lesssim 2\sigma$) for
the same.

##  I Introduction

To date, the majority of cosmological observations Adam  _et al._ (2016); Ade
_et al._ (2016a) related to the space-time of the Universe can be explained by
General Relativity (GR), particularly in the smaller scales. To explain the
observed late time acceleration Riess  _et al._ (1998); Perlmutter  _et al._
(1999) with metric based GR framework, one can: i) introduce a cosmological
constant Weinberg (1989); Bull  _et al._ (2016), ii) introduce dark energy
(DE) such as a dynamical scalar field (e.g., quintessence Caldwell  _et al._
(1998), k essence Chiba  _et al._ (2000); Armendariz-Picon  _et al._ (2001))
or (iii) modify the gravitational coupling differing from GR at cosmological
distance, known as the modified gravity (MG) theory Tsujikawa (2010); Nojiri
and Odintsov (2011); Clifton  _et al._ (2012); Joyce  _et al._ (2016).

To test a large number of MG models, a generalized platform is necessary. The
degenerate higher-order scalar-tensor (DHOST) theory is claimed to be the most
general class of a scalar-tensor theory, where a propagating scalar degree of
freedom (DOF) is added on the top of two tensor DOF of metric based GR given
under the general-covariance (see Appendix A for more details) Langlois and
Noui (2016); Crisostomi  _et al._ (2016); Ben Achour  _et al._ (2016a);
Motohashi  _et al._ (2016); Ben Achour  _et al._ (2016b) Kobayashi (2019). The
DHOST theory includes many popular modified gravity models such as Brans-Dicke
theory Brans and Dicke (1961), $f(R)$ gravity De Felice and Tsujikawa (2010);
Nojiri and Odintsov (2011), covariant Galileon Deffayet  _et al._ (2009);
Neveu  _et al._ (2013, 2017), Horndeski Horndeski (1974); Kobayashi  _et al._
(2011), transforming gravity Zumalacárregui and García-Bellido (2014), and
GLPV theory Gleyzes  _et al._ (2015). A large class of DHOST theories produces
gravitational waves with velocity ($c_{g}$) equal to the velocity of light
($c$) in agreement with recent observations. These are referred to as viable
DHOST$c_{g}^{2}=c^{2}$ theories in this article.

In order to explain the small scale observations, the additional degrees of
freedom (DOF) should be screened through a non-linear mechanism, and GR should
be recovered on the small scale, known as the so-called screening mechanism.
The Vainshtein screening Vainshtein (1972) is useful for the higher-order MG
theories Babichev and Deffayet (2013); Brax  _et al._ (2004). In Vainshtein
screening, the gravitational potentials for the DHOST$c_{g}=c$ are modified
inside the matter sources Crisostomi and Koyama (2018); Langlois  _et al._
(2018a); Bartolo  _et al._ (2018); Dima and Vernizzi (2018); Hirano  _et al._
(2019); Crisostomi  _et al._ (2019). Horndeski theory, a subset of DHOST
theory, recovers GR in the small scale approximation Kimura  _et al._ (2012);
Narikawa  _et al._ (2013); Koyama  _et al._ (2013); Kase and Tsujikawa (2013).
However, Vainshtein screening for GLPV theory, another subset of DHOST theory,
breaks down inside the matter Kobayashi  _et al._ (2015); Koyama and Sakstein
(2015); Saito  _et al._ (2015); Sakstein (2015a, b); Jain  _et al._ (2016);
Babichev  _et al._ (2016). DHOST theory has not been studied extensively on
the galaxy cluster scale, other than a few bounds on $\Xi_{1}$ Sakstein  _et
al._ (2016); Salzano  _et al._ (2017). However, as shown in Saito  _et al._
(2015), one must always have $\Xi_{1}>-1/6$ to guarantee a stable static
solution for stars, and the upper bound $\Xi_{1}<7\times 10^{-3}$ to fulfill
the consistency of the minimum mass for hydrogen burning in stars with the
lowest mass red dwarf Sakstein (2015a, b). Which will be very useful in order
to compare with the constraints that can be obtained from the galaxy cluster
scales.

Given the task at hand, we utilize the XMM-Newton Cluster Outskirts Project
(X-COP) data products (Eckert  _et al._ , 2017; Ettori  _et al._ , 2019;
Eckert  _et al._ , 2019; Ghirardini  _et al._ , 2019), which consists of 12
clusters with well-observed X-ray emission and high signal to noise ratio in
the Planck Sunyaev-Zel’dovich (SZ) survey (Ade  _et al._ , 2016b), essentially
providing both ICM temperature and pressure data over the large radial range
of $0.2\,{\rm Mpc}\lesssim r\lesssim 2\,{\rm Mpc}$. Such a compilation is in
fact very helpful to assess the mass profiles out to a large radial range and,
as in the current work, subsequently test modifications to the gravitational
potential.

A similar approach relying on the same galaxy cluster phenomenology has been
utilized in several previous analyses (Sakstein  _et al._ , 2016) (see in
particular (Terukina  _et al._ , 2014; Wilcox  _et al._ , 2015) for tests of
Chameleon screening). In (Sakstein  _et al._ , 2016) utilize both X-ray and
weak lensing profiles of 58 stacked clusters, obtaining constraints on the
modifications to gravitational potential which point towards consistency with
GR. On the contrary, in this work for the very first time (as far as we are
aware) we exploit individual clusters with well-observed ICM to finally obtain
a joint constraint, while simultaneously assessing a possible redshift
dependence. More recently, (Pizzuti  _et al._ , 2020) have presented
generalized frameworks for testing these scenarios, also including the beyond
Horndeski and DHOST scenarios, essentially constraining the Vainshtein
screening.

The paper is organized as follows. We provide a brief description of the DHOST
theory in Section II, while the modeling of the hydrostatic equilibrium in the
case of a cluster of galaxies is summarized in Section III. The datasets
utilized and method is described in Section IV. Finally we present and discuss
our results in Section VI. In the current analysis, we assume $H_{0}=70$ km/s
Mpc-1and $\Omega_{\rm m}=0.3$, when computing the critical density $\rho_{\rm
c}(z)=3H^{2}(z)/8\pi G_{\rm N}^{\rm eff}$111For the case of GR $G_{\rm N}^{\rm
eff}\to G_{\rm N}$, please see Section II for details., where
$H^{2}(z)/H^{2}_{0}=\Omega_{\rm m}(1+z)^{3}+(1-\Omega_{\rm m})$.

##  II Perturbed gravitational forces of DHOST in small spatial scale

Gravity is well tested at small scales. We expect to recover GR within the
Vainstein radius (including inside the object), and full extended theory
outside. However, for the DHOST theory, GR is not recovered everywhere inside
the Vainstein radius, i.e., Vainstein screening breaks down inside the object.
The gravitational potentials for the DHOST theory in the galaxy scales
modifies to Langlois  _et al._ (2018b); Crisostomi and Koyama (2018); Bartolo
_et al._ (2018); Dima and Vernizzi (2018)

$\displaystyle\frac{d\Phi(r)}{dr}$ $\displaystyle=$ $\displaystyle\frac{G_{\rm
N}^{\rm eff}{\cal M}(r)}{r^{2}}+\Xi_{1}G_{\rm N}^{\rm eff}{\cal
M}^{\prime\prime}(r)\,,$ (1) $\displaystyle\frac{d\Psi(r)}{dr}$
$\displaystyle=$ $\displaystyle\frac{G_{\rm N}^{\rm eff}{\cal
M}(r)}{r^{2}}+\Xi_{2}\frac{G_{\rm N}^{\rm eff}{\cal
M}^{\prime}(r)}{r}+\Xi_{3}\,G_{\rm N}^{\rm eff}{\cal M}^{\prime\prime}(r)\,,$
(2)

with the modified or effective Newton’s constant $G_{\rm N}^{\rm eff}$ defined
by the expression

$\displaystyle\frac{1}{8\pi G_{\rm N}^{\rm eff}}$ $\displaystyle=$
$\displaystyle{2F-2{F_{X}}X-\frac{3}{2}A_{3}X^{2}}\equiv 2F(1+\Xi_{0})\,,$ (3)

and the dimensionless coefficients

$\left.\begin{aligned} \Xi_{1}&=-\frac{(4{F_{X}}-XA_{3})^{2}}{16FA_{3}}\,,\\\
\Xi_{2}&=-\frac{2{XF_{X}}}{F}\,,\\\
\Xi_{3}&=\frac{16{F_{X}}^{2}-A_{3}^{2}X^{2}}{16A_{3}F}\,,\\\
\Xi_{0}&=-\frac{F_{X}X}{F}-\frac{3}{4}\frac{A_{3}X^{2}}{F}\,.\end{aligned}\right.$
(4)

From (4), one can derive a consistency relation,
$\Xi_{2}=(2/\Xi_{1})\left(\Xi_{3}^{2}-\Xi_{1}^{2}\right)$, which allows to
reduce to $(\Xi_{1},\Xi_{3})$ the number of parameters regulating the
deviations from $\Lambda$CDM in the perturbations regime Langlois  _et al._
(2018a).

One can recover GR by setting $F=1/2\kappa$ with $\kappa=8\pi G_{\rm
N}c^{-4}$, and $P,Q,A_{3}=0$ (see Appendix A), which would lead to
$\Xi_{0,1,2,3}=0$.

In terms of the cosmological effective parameters (EFT parameters) Langlois
_et al._ (2017),

$\left.\begin{aligned}
\Xi_{1}&=-\frac{(\alpha_{\text{H}}+\beta_{1})^{2}}{2(\alpha_{\text{H}}+2\beta_{1})}\,,\quad\Xi_{2}=\alpha_{\text{H}}\,,\quad\\\
\Xi_{3}&=-\frac{\beta_{1}(\alpha_{\text{H}}+\beta_{1})}{2(\alpha_{\text{H}}+2\beta_{1})}\,,\quad\Xi_{0}=-\alpha_{\text{H}}-3\beta_{1}\,.\end{aligned}\right.$
(5)

Therefore, $G_{\rm N}^{\rm eff}\equiv G_{\rm N}^{\rm eff}(\alpha_{\rm
H},\beta_{1})$ and $\Xi_{\mu}\equiv\Xi_{\mu}(\alpha_{\rm H},\beta_{1})$ and
can be explicitly written as,

$\left.\begin{aligned} G_{\rm N}^{\rm eff}&=[16\pi F(1+\Xi_{0})]^{-1}\\\
&=[8\pi M_{\star}^{2}(1+\Xi_{0})]^{-1}\\\ &={\gamma}_{\rm N}G_{\rm
N}/(1-\alpha_{\rm H}-3\beta_{1})\,,\end{aligned}\right.$ (6)

where the effective Planck mass $M_{\star}=M_{\rm Pl}/{\gamma}_{\rm N}$.
Hereafter, we write $G_{\rm N}^{\rm eff}=\tilde{\gamma}_{\rm N}G_{\rm N}$ for
brevity. Note that it is always possible to convert the parameterization from
the effective ones to the physical ones after the analysis is performed. We
elaborate on this in Appendix B.

As the current modification is a linear order perturbation to the metric, the
EFT parameters depend only on the background. Thereby, the $G_{\rm N}^{\rm
eff}(t)$ and $\Xi_{i}(t)$ are solely functions of time.

##  III ICM Pressure profile of galaxy cluster

As per the standard model, the galaxy clusters consists of Dark Matter (DM)
($\sim$ 85% of the total mass), the hot ionized hydrogen and helium gas called
intracluster medium (ICM) ($\sim$12%), and visible stars and galaxies
($\sim$3%). The galaxy clusters are expected to retain all the ICM and stars
accreted since the formation epoch in the deep gravitational well created by
DM. On one hand, the hot ICM is visible in the X-ray band through thermal
bremsstrahlung and line emission. On the other hand, the ICM can be detected
at millimetre wavelengths through the distortion of the cosmic microwave
background (CMB) induced by inverse Compton scattering, called the SZ effect.
Since ICM traces the main baryonic component, X-ray and SZ observations are
crucial for testing gravity in the clusters and filaments, which is however a
challenging task.

Among the possible approaches to infer the total mass, one possible method is
to assume that the ionized gas is in hydrostatic equilibrium with the
gravitational potential, $\Phi$, (mostly created by DM (Kravtsov and Borgani,
2012)). In the hydrostatic equilibrium,

$\frac{1}{\rho_{\rm gas}(r)}\frac{{\rm d}P_{\rm gas}(r)}{{\rm d}r}=-\frac{{\rm
d}\Phi(r)}{{\rm d}r}\,.$ (7)

where we have implicitly assumed spherical symmetry. $P_{\rm gas}(r)$ is the
gas pressure radial profile while the gas density is:

$\displaystyle\rho_{\rm gas}(r)$ $\displaystyle=$ $\displaystyle n_{\rm
gas}(r)\,\mu m_{p}\,,$ (8) $\displaystyle\sim$ $\displaystyle 1.8n_{e}(r)\,\mu
m_{p}\,,$ (9)

where $n_{\rm gas}(r)$ is the gas density, the sum of the electron and proton
number densities, i.e., $n_{\rm gas}(r)=n_{e}(r)+n_{p}(r)$, $\mu$ is the mean
molecular weight in a.m.u. ($\mu=0.61$ is the corresponding mean molecular
weight in atomic mass unit), and $m_{p}$ is the proton mass ($\sim$ 1 atomic
mass unit, $m_{u}$) therefore $n_{p}\approx 0.8n_{e}$ which gives $n_{\rm
gas}(r)\approx 1.826\,n_{e}(r)$.

Among the possible models for $n_{e}(r)$, such as the $\beta$ profile, or
double $\beta$ profile, or the more recent Vikhlinin profile, we choose the
latter Vikhlinin  _et al._ (2006). We use the simplified Vikhlinin222Here we
have utilized the Vikhlinin profile with only 6 parameters, neglecting a
second part of the parametric model which adds 3 more parameters. We comment
on the implications for this at in Appendix C. parametric model setting,

$\frac{n_{e}(r)}{n_{0}}=\frac{(r/r_{c})^{-\alpha/2}[1+(r/r_{s})^{\gamma}]^{-\varepsilon/(2\gamma)}}{[1+(r/r_{c})^{2}]^{(3/2)\beta-\alpha/4}}$
(10)

where $n_{0}$ is the cluster central density, $r_{s}$ is the transition radius
at which the logarithmic slope changes, and $r_{c}$ is the core radius. The
$\beta$ and $\epsilon$ parameters give the inner and outer profile slope,
respectively. The $\gamma$ parameter gives the width of the transition in the
profile. We fix $\gamma=3$ as suggested in Vikhlinin  _et al._ (2006);
Vikhlinin (2006), so the electron density parameter space is given by
${\Theta_{e}}=\\{n_{0},\alpha,\beta,\varepsilon,r_{\rm c},r_{\rm s}\\}$.
Moreover, we limit $\epsilon<5$ to avoid nonphysical sharp features in the
electron density profile, which we however find to not make a major difference
on the mass estimates once rest of the parameters are marginalized upon.

The right side of Equation 7 is governed by the gravitational potential, which
depends on the assumed theory of gravitation. In the standard GR scenario we
have,

$\displaystyle\frac{d\Phi(r)}{dr}$ $\displaystyle=$ $\displaystyle\frac{G_{\rm
N}}{r^{2}}{M}_{\rm HSE}(r).$ (11)

For the DHOST theory elaborated in Equation 1, the gravitational potential can
be written as,

$\frac{d\Phi(r)}{dr}=\frac{G_{\rm N}^{\rm eff}M_{\rm
HSE}(r)}{r^{2}}+\Xi_{1}G_{\rm N}^{\rm eff}\frac{d^{2}M_{\rm HSE}(r)}{dr^{2}}.$
(12)

In Equation 12, $G_{\rm N}^{\rm eff}$ is the effective gravitational constant
which is related to the Newton one $G_{\rm N}$ in the DHOST model through
Equation 6. One could ideally write $G_{\rm N}^{\rm eff}=\tilde{\gamma}_{\rm
N}G_{\rm N}$, where $\tilde{\gamma}_{\rm N}$ is a redshift dependent function.
However, it is evident that $\tilde{\gamma}_{\rm N}$ can not be constrained by
the data on the pressure profile alone since it only enters as a
multiplicative term fully degenerate with the total mass. Also, for viable
DHOST theories, it is expected that $G_{\rm N}^{\rm eff}$ does not deviate
from $G_{\rm N}$ so we perform the analysis fixing $\tilde{\gamma}_{\rm N}=1$,
and then $\tilde{\gamma}_{\rm N}\neq 1$.

Under the assumption that the cluster mass is dominated by the dark matter
component, we model the mass density using the NFW profile Navarro  _et al._
(1996) given as,

$\rho(r)=\frac{\rho_{s}}{(r/r_{s})(1+r/r_{s})^{2}}$ (13)

where $\rho_{s}$ is a characteristic density, and $r_{s}$ the radius where the
logarithmic slope $s=d\ln{\rho}/d\ln{r}$ takes the isothermal value $s=-2$.
Now the mass profile can be straightforwardly obtained and it is conveniently
rewritten as

$M(<r)=M_{500}\frac{\ln{(1+c_{500}x)}-c_{500}x/(1+c_{500}x)}{\ln{(1+c_{500})}-c_{500}/(1+c_{500})}$
(14)

with $x=r/R_{500}$, $c_{500}=R_{500}/r_{s}$333$R_{500}$ in the DHOST case
already takes into account the variation due to modification to gravity on the
background through $\rho_{\rm c}(z)$, hence having
$R_{500}^{3}\sim\tilde{\gamma}_{\rm N}M_{500}$. the halo concentration, and

$M_{500}=500\,\frac{4}{3}\pi\rho_{\rm c}(z)R_{500}^{3}\ .$ (15)

From Equation 14, we obtain,

$M^{\prime\prime}(r)=\frac{M_{500}}{R_{500}^{2}}\frac{c_{500}^{2}(1-c_{500}x)(1+c_{500}x)^{-3}}{\ln{(1+c_{500})}-c_{500}/(1+c_{500})}\
,$ (16)

which can be plugged into Equation 12 together with Equation 14 to get the
modified gravitational potential. Now, the parameters for the mass profile are
$\Theta_{\rm M}=\\{M_{500},c_{500}\\}$. One can get the SZ pressure profile
for the galaxy cluster by integrating (7) and using (1),

$P^{\rm th}(r)=P^{\rm th}(0)-1.8\mu m_{\rm p}\int_{0}^{r}n_{\rm
e}(\tilde{r})\left[\frac{G_{\rm N}^{\rm eff}{M}_{\rm
HSE}(\tilde{r})}{\tilde{r}^{2}}+\Xi_{1}G_{\rm N}^{\rm
eff}{M^{\prime\prime}}_{\rm HSE}(\tilde{r})\right]{\rm d}\tilde{r}$ (17)

where, $G_{\rm N}^{\rm eff}=\tilde{\gamma}_{\rm N}\times G_{\rm N}$ and
$M_{\rm HSE}$ are modeled using the NFW profile Equation 14. Note that here
$P^{\rm th}(0)$ is not a free parameter and can be reconstructed form the
integration assuming $P^{\rm th}(r\to\infty)\sim 0$. Effectively, the DHOST
modification to the standard GR case, is dictated by the combination of two
parameters $\\{\Xi_{1},\tilde{\gamma}_{\rm N}\\}$, which we hereafter denote
as $\Theta_{\rm DHOST}$. The formalism adopted here is in fact termed as
Backward method, where a given parametric profile is assumed for the mass
model and the pressure is obtained by the integration of the same. In
contrast, assuming a parametric profile for the pressure is termed as Forward
method, in which case assessing the modifications to gravitational potential
is not possible.

##  IV Data and Likelihood

We utilize the collection of clusters within the XMM-Newton Cluster Outskirts
Project (X-COP)444The datasets are publicly available at the following link:
https://dominiqueeckert.wixsite.com/xcop/about-x-cop catalog (Eckert  _et al._
, 2017) with joint X-ray and millimeter observations. The compilation consists
of 12 massive ($10^{14}\,M_{\odot}\lesssim M_{500}\lesssim
10^{15}\,M_{\odot}$) galaxy clusters in the redshift range $0.04<z<0.1$,
selected for high signal-to-noise ratios (${\rm S/N}>12$) in the Planck SZ
survey Ade  _et al._ (2014). The physical observables of interest are, $i)$
Electron density, $ii)$ Temperature of gas, in the X-ray and the $iii)$
Pressure of the gas observed in the SZ. These observables indeed present a
suitable scenario to test the formalism described in Section III and the
modifications of gravity on cluster scales.

In Figure 1, we present the collection of the 12 clusters data currently
analyzed. Firstly, we re-scale the data for the self-similar normalization
(i.e, the data were provided as a function of $R/R_{500}$, for each cluster),
which they were provided with and show the radial profiles.

Figure 1: Top: Electron density radial profiles obtained using the L1
regularization method (Ghirardini  _et al._ , 2019). Middle: Temperature of
the X-ray emitting gas. Bottom: The SZ electron pressure profile. This figure
is comparable to the Fig.1 of Ettori  _et al._ (2019), except that we have
removed the self-similar scaling, as reported therein. Also note that they
show the reconstructed electron density profiles.

In order to constrain the DHOST modification parameters ($\Theta_{\rm
DHOST}$), the halo mass profile parameters ($\Theta_{M}$), and the electron
density parameters (${\Theta}_{e}$), we perform a joint fit to the measured
$n_{e}(r)$, $P_{\rm X}(r)$ and $P_{\rm SZ}(r)$ data from the X-COP sample. The
total number of parameters is then,

${\Theta}=\Theta_{e}\cup{\Theta}_{\rm M}\cup\Theta_{\rm DHOST}.$ (18)

accounting for a total of 10 parameters in the DHOST scenarios and 8 within GR
case, having $\Theta_{\rm DHOST}\equiv\\{0.0,1.0\\}$. Similar to Ettori  _et
al._ (2019); Ghirardini  _et al._ (2018), we define the likelihood function
for the observed data denoted as ($\rm obs$) and for the model ($\rm mod$) as
in the following way,

$\displaystyle-2\ln{{\cal{L}}}$ $\displaystyle=$ $\displaystyle({\bf P}^{\rm
obs}_{\rm SZ}-{\bf P}^{\rm mod}_{\rm SZ})\Sigma_{\rm TOT}^{-1}({\bf P}^{\rm
obs}_{\rm SZ}-{\bf P}^{\rm mod}_{\rm SZ})^{T}{+\ln\mid\Sigma_{\rm TOT}\mid}$
(19) $\displaystyle+$ $\displaystyle\sum_{j=1}^{N_{P_{\rm
X}}}\left[\frac{(P_{{\rm X},i}^{\rm obs}-P_{{\rm X},i}^{\rm
mod})^{2}}{\sigma_{P_{{\rm X},i}}^{2}+\sigma_{P,{\rm
int}}^{2}}+\ln({\sigma_{P_{{\rm X},i}}^{2}+\sigma_{P,{\rm int}}^{2}})\right]$
$\displaystyle+$ $\displaystyle\sum_{i=1}^{N_{\rm
n_{e}}}\left[\frac{(n_{e,i}^{\rm obs}-n_{e,i}^{\rm
mod})^{2}}{\sigma_{n_{e,i}}^{2}}\right]$

where the first term accounts for the co-varying SZ pressure data, with the
covariance matrix $\Sigma_{\rm TOT}=\Sigma_{\rm P}+\Sigma_{{\rm
P,int}}$555$\Sigma_{{\rm P,int}}$ is a diagonal matrix of $\sigma_{P,{\rm
int}}^{2}$, accounting for an additional intrinsic scatter for the $P_{\rm X}$
and $P_{\rm SZ}$ data points., $\Sigma_{\rm P}$ being the covariance matrix of
the $P_{\rm SZ}$ data itself. The second and the third term account for the
X-ray pressure and the electron density data, respectively. Note that
following (Ghirardini  _et al._ , 2018; Ettori  _et al._ , 2019), we include
an additional intrinsic scatter to the measurement uncertainty of the pressure
and the temperature datasets. Our Likelihood is practically identical to the
one implemented in Ettori  _et al._ (2019), except that we have utilized the
publicly available electron density data instead of the emissivity
$(\epsilon)$ data utilized therein. The electron density data taken into
account in our analysis is obtained through the L1 regularization method
(Ghirardini  _et al._ , 2019), in contrast to the smooth reconstructed
electron density profiles utilized in Ettori  _et al._ (2019). Their
reconstructed profiles were obtained through the multiscale fitting Eckert
_et al._ (2016) procedure to the emissivity data. The multiscale fitting of
the emissivity data essentially provides smooth reconstructed electron density
with minimal features and in this way the reconstructed electron density
profiles show limited features resulting in a better estimate of the
hydrostatic mass, having better control over the gradient. The L1
regularization data we utilize here shows more pronounced features666The gas
clumping in the electron density data should be accounted for while obtaining
the smoothed reconstructed profiles. However, the possible bias due to the
clumping is mild Eckert  _et al._ (2015).  and represents the electron density
observations better to study the shape of the profile (Ghirardini  _et al._ ,
2019). However, we proceed to utilize the data here, expecting a mild to
moderate deviation from the hydrostatic masses reported in (Ettori  _et al._ ,
2019). This is partly due to the fact that, in (Ghirardini  _et al._ , 2019),
it has been shown that the two methods agree incredibly well and show an
utmost scatter of the order $\sim 5\%$ at a given radius, at least for the
clusters at hand in the current compilation. And that using the simplified
Vikhlinin profile (Equation 10) provides smooth reconstructions for the
electron density and varies even more minimally with respect to the multiscale
fitting reconstruction. However, the L1 based electron density data also have
mild correlations among the data points777 Dominique Eckert and Stefano Ettori
in private communication., which we have not taken into account here and could
potentially have an impact on the final mass estimates. In any case, as the
current analysis is performed to assess the deviations from GR to the DHOST
scenarios, we expect the use of $n_{e}$ data as here, instead of the
$\epsilon$ data, to equivalently affect the GR and the DHOST scenarios. This
we indeed verify a posteriori, finding minimal correlations between the
$\Theta_{\rm M}$ and $\Theta_{e}$ parameters.

We also utilize the likelihood using the $T_{\rm X}$ data to perform
validations, which is the direct observable while $P_{\rm X}$ is obtained by
$T_{\rm X}$ data and the measured electron density. In this case, the second
term in the Equation (19) is replaced with the $T_{\rm X}$ data and the rest
remains unchanged. As shown in (Ghirardini  _et al._ , 2018), a small
intrinsic scatter ($\sigma_{\rm int}$) is modeled up on $\log(P)$ such that
the intrinsic error on the pressure can be written as $\sigma_{P,{\rm
int}}=P\,\sinh({\sigma_{\rm int}})$. In our analysis we utilize the
$\sinh({\sigma_{\rm int}})$ as a free parameter888It is usually convenient to
sample on the $\sinh({\sigma_{\rm int}})$ as a free parameter, to avoid loss
of probability for the posterior of parameter $\sigma_{\rm
int}\xrightarrow{}0$ (Hogg and Foreman-Mackey, 2018).  and the intrinsic
scatter is propagated to the temperature profile as $\sigma_{T,{\rm
int}}=T\,\sinh({\sigma_{\rm int}})$. Therefore, the intrinsic scatter is
included to the parameter array $\Theta$ in Equation 18, as an additional
parameter.

Finally, the posterior defined in Equation (19), is utilized to perform a
Bayesian analysis through MCMC sampling. We use the
emcee999http://dfm.io/emcee/current/ package (Foreman-Mackey  _et al._ , 2013;
Hogg and Foreman-Mackey, 2018), which implements an affine-invariant ensemble
sampler. We utilise publicly available
ChainConsumer101010https://github.com/Samreay/ChainConsumer/tree/Final-Paper.
package (Hinton, 2016), to perform analysis of the chains and plot the
contours. To asses the evidence for the DHOST modification we compute the
Bayesian evidence Trotta (2008, 2017); Heavens  _et al._ (2017a), through the
MCEvidence package Heavens  _et al._ (2017b)111111We utilized the MCEvidence
package publicly available at https://github.com/yabebalFantaye/MCEvidence..
We impose uniform flat priors on all the parameters, in particular we impose
$\Xi_{1}\in\\{-2.0,2.0\\}$ and $\tilde{\gamma}_{\rm N}\in\\{0.001,2.0\\}$, for
the DHOST parameters. We also separately test the effects of different
initialization of the mass profile parameters for the walkers in the affine-
invariant procedure, finding no implications for the posteriors. Note that the
information on the prior volume is important in assessing the Bayesian
evidence. Through the Bayes’ rule the posterior distribution w.r.t. the
parameters $(\Theta)$ of the given model $\mathcal{M}(\Theta)$ and the
observations $\mathcal{D}$, can be written as,

$p(\Theta|\mathcal{D},\mathcal{M})=\dfrac{p(\mathcal{D}|\Theta,\mathcal{M})\pi(\Theta|\mathcal{M})}{p(\mathcal{D}|\mathcal{M})},$
(20)

where $\pi(\Theta|\mathcal{M})$ is the prior volume and
$p(\mathcal{D}|\mathcal{M})$ is the Bayesian ‘evidence’ $(\cal{B})$. A
comparison of the evidence can be utilized to assess the preference of a given
model $\mathcal{M}_{1}(\Theta_{1})$ over the base model, which in our case are
the DHOST and GR, respectively. As is the usual practice we contrast the
Bayesian evidence using the Jeffrey’s scale (Jeffreys, 1961), where
$\Delta\log(\cal{B})$ $<2.5$ and $\Delta\log(\cal{B})$ $>5$, imply a weak and
strong preference for the extended model against the base model, respectively.

##  V Results and Discussion

Figure 2: Probability density for the distribution of the parameter $\Xi_{1}$.
We perform a simple Gaussian kernel density estimation to obtain the smooth
profiles, which are over-plotted on the distributions. The vertical dashed
line marks the GR ($\Xi_{1}=0$) case.

We begin by performing the analysis for the standard GR case, aiming to
reproduce the results of E19. The constraints obtained in our analysis are
shown in the first three columns of Table 1. We find a good agreement with the
constraints reported in E19, however, note that here we have reported
asymmetric uncertainties accounting for the maximum likelihood, in contrast to
the symmetric dispersion between the $84^{\rm th}$ and $16^{\rm th}$
percentiles presented therein. In general, we find a good agreement within the
order of $\lesssim 1.5\sigma$, for the mass parameter $M_{500}$, also with
very similar uncertainty estimates. We stress that a difference of the order
$\sim 1.5\sigma$ is yet within the usual expectations of a bias using
different methods or datasets of estimating $M_{500}$. The difference between
the analysis in E19 and here is within the modeling of the electron density
data, where they utilize the smooth reconstructed profiles from the multiscale
fitting, whereas we utilize the publicly available L1 regularization based
data points (see (Ameglio  _et al._ , 2009; Croston  _et al._ , 2006)). The
presence of clumpiness, if assumed to be smooth essentially yields an
underestimate of the total mass (Mathiesen  _et al._ , 1999; Ettori  _et al._
, 2013), and should indicate a systematic bias in our mass estimates. However,
as can be seen in Figure 8, our mass estimates do not show such a systematic
behavior allowing us to validate the formalism implemented here. Therefore,
owing to the agreement with the analysis in E19 and no systematic behavior in
mass constraints we find our results to be consistent with those presented in
E19 and proceed with further analysis of comparing the GR and modified gravity
scenarios, where any difference from the earlier results will be equivalently
present in both cases and is not expected to bias our inferences for the DHOST
scenario. We present a more detailed comparison of masses both in the GR and
the DHOST cases in Appendix C.

Table 1: Constraints obtained for the parameters of the mass profiles in both the GR and the DHOST scenarios. We show the maximum likelihood statistics as the $68\%$ C.L. limits. Here we have utilized the $P_{\rm SZ}+P_{\rm X}$ datasets. In the last column, we show the Bayesian evidence for the DHOST modification w.r.t GR, where a negative value indicates that GR is preferred. The four non-NFW clusters are denoted with the superscript † and additionally, the cluster A2319 carries ‡, representing the large non-thermal effects reported in (Eckert  _et al._ , 2019).  Cluster | GR | DHOST |
---|---|---|---
| $c_{500}$ | $M_{500}\,$ | $R_{500}\,$ | $c_{500}$ | $M_{500}\,$ | $R_{500}\,$ | $\Xi_{1}$ | $\tilde{\gamma}_{\rm N}\times M_{500}\,$ | $\Delta\log(\cal{B})$
|  | $[10^{14}\,M_{\odot}]$ | $[{\rm Mpc}]$ |  | $[10^{14}\,M_{\odot}]$ | $[{\rm Mpc}]$ |  | $[10^{14}\,M_{\odot}]$ |
| A85
---
(z=0.0555)
$2.05^{+0.09}_{-0.06}$ | $6.14^{+0.14}_{-0.21}$ | $1.270^{+0.010}_{-0.015}$ | $1.16^{+0.57}_{-0.32}$ | $4.10^{+3.43}_{-0.63}$ | $1.292^{+0.017}_{-0.030}$ | $0.30^{+0.11}_{-0.27}$ | $6.46^{+0.26}_{-0.45}$ | $-3.7$
| A644†
---
(z=0.0704)
$4.22^{+0.31}_{-0.17}$ | $4.93^{+0.25}_{-0.18}$ | $1.175^{+0.020}_{-0.015}$ | $6.76^{+0.28}_{-0.33}$ | $3.80^{+5.13}_{-0.71}$ | $0.980^{+0.028}_{-0.030}$ | $-1.04^{+0.18}_{-0.19}$ | $2.85^{+0.24}_{-0.26}$ | $46.5$
| A1644†
---
(z=0.0473)
$0.95^{+0.10}_{-0.10}$ | $3.00^{+0.17}_{-0.15}$ | $1.003^{+0.019}_{-0.017}$ | $2.89^{+0.23}_{-0.20}$ | $2.03^{+3.77}_{-0.67}$ | $0.844^{+0.020}_{-0.027}$ | $-0.837^{+0.119}_{-0.090}$ | $1.78^{+0.13}_{-0.16}$ | $39.9$
| A1759
---
(z=0.0622)
$3.08^{+0.15}_{-0.10}$ | $4.59^{+0.18}_{-0.12}$ | $1.150^{+0.015}_{-0.010}$ | $3.71^{+0.47}_{-0.35}$ | $2.88^{+2.69}_{-0.81}$ | $1.101^{+0.032}_{-0.035}$ | $-0.169^{+0.111}_{-0.090}$ | $4.03^{+0.34}_{-0.39}$ | $-2.9$
| A2029
---
(z=0.0773)
$3.14^{+0.12}_{-0.17}$ | $7.84^{+0.33}_{-0.26}$ | $1.369^{+0.019}_{-0.015}$ | $3.31^{+0.49}_{-0.84}$ | $6.0^{+4.9}_{-1.6}$ | $1.352^{+0.089}_{-0.016}$ | $-0.04^{+0.19}_{-0.12}$ | $8.04^{+0.96}_{-0.88}$ | $-3.0$
| A2142
---
(z=0.0909)
$2.25^{+0.10}_{-0.12}$ | $8.30^{+0.33}_{-0.28}$ | $1.389^{+0.017}_{-0.017}$ | $2.86^{+0.33}_{-0.29}$ | $4.55^{+6.06}_{-0.92}$ | $1.326^{+0.040}_{-0.024}$ | $-0.203^{+0.101}_{-0.079}$ | $7.21^{+0.65}_{-0.40}$ | $0.4$
| A2255†
---
(z=0.0809)
$0.68^{+0.13}_{-0.10}$ | $5.02^{+0.31}_{-0.26}$ | $1.180^{+0.023}_{-0.021}$ | $2.44^{+0.15}_{-0.21}$ | $5.2^{+2.4}_{-2.1}$ | $0.953^{+0.046}_{-0.043}$ | $-1.1^{+0.26}_{-0.32}$ | $2.66^{+0.38}_{-0.36}$ | $28.2$
| A2319${}^{\dagger,\,\ddagger}$
---
(z=0.0557)
$3.40^{+0.13}_{-0.09}$ | $7.15^{+0.16}_{-0.09}$ | $1.336^{+0.010}_{-0.006}$ | $5.14^{+0.13}_{-0.16}$ | $8.83^{+3.63}_{-2.26}$ | $1.151^{+0.020}_{-0.016}$ | $-0.827^{+0.108}_{-0.076}$ | $4.57^{+0.23}_{-0.19}$ | $109.6$
| A3158
---
(z=0.0597)
$1.81\pm 0.12$ | $4.21^{+0.19}_{-0.14}$ | $1.119^{+0.016}_{-0.012}$ | $2.62^{+0.38}_{-0.56}$ | $2.46^{+2.90}_{-0.62}$ | $1.054^{+0.057}_{-0.029}$ | $-0.23^{+0.15}_{-0.18}$ | $3.51^{+0.57}_{-0.31}$ | $-0.4$
| A3266
---
(z=0.0589)
$0.93^{+0.10}_{-0.10}$ | $9.90^{+0.57}_{-0.59}$ | $1.489^{+0.027}_{-0.030}$ | $0.71^{+0.17}_{-0.20}$ | $6.3^{+4.0}_{-1.8}$ | $1.455^{+0.045}_{-0.055}$ | $0.100^{+0.137}_{-0.079}$ | $9.23^{+0.87}_{-0.99}$ | $-1.3$
| RXC1825
---
(z=0.0650)
$2.16^{+0.15}_{-0.12}$ | $4.11^{+0.15}_{-0.13}$ | $1.108^{+0.013}_{-0.012}$ | $1.54^{+0.44}_{-0.43}$ | $2.68^{+2.49}_{-0.46}$ | $1.130^{+0.016}_{-0.018}$ | $0.17^{+0.17}_{-0.13}$ | $4.37^{+0.17}_{-0.22}$ | $-0.3$
| ZW1215
---
(z=0.0766)
$1.32^{+0.11}_{-0.14}$ | $7.82^{+0.51}_{-0.50}$ | $1.368^{+0.029}_{-0.029}$ | $1.98^{+0.51}_{-0.84}$ | $7.23^{+4.31}_{-2.85}$ | $1.331^{+0.041}_{-0.076}$ | $-0.21^{+0.27}_{-0.18}$ | $7.14^{+0.65}_{-1.18}$ | $-4.5$

###  V.1 Constraints on the DHOST modification

The pressure profile for the DHOST gravity is estimated following Equation 17,
which now accounts for the modified background evolution and the small scale
clustering effects through the effective parameters $\tilde{\gamma}_{\rm N}$
and $\Xi_{1}$, respectively. We show the constraints on the mass and the DHOST
parameters in Table 1, under the column titled DHOST. To begin with, we notice
that the clusters A644, A1644, A2319 and A2255 show large negative values for
$\Xi_{1}$ parameter, significantly different from the GR ($\Xi_{1}=0$)
expectation. However, it has already been shown in E19 that these same four
clusters present a mild preference for mass profiles other than the NFW one,
with a mild-to-moderate evidence for either the King approximation to
Isothermal sphere (A2225, A2319), the Burkert (A644), or the Hernquist (A1644)
profiles. Moreover, these 4 clusters are also outliers with respect to the
$c_{200}-M_{200}$ scaling relation. Therefore, we do not immediately take our
constraints for these clusters as evidence for modification of gravity, but
warn the reader that any inferences should be made with caution. Hereafter we
refer to these four clusters as ‘non-NFW’ clusters for the ease of discussion
and represent them with a superscript † in the tables. This clearly indicates
a degenerate scenario between an assumption of the mass profile in the GR case
and a modification to the gravity itself, such as the DHOST theory in the
current work.

Figure 3: Constraints on $\Xi_{1}$ and the ratios of $\tilde{\gamma}_{\rm
N}\times M_{500}^{\rm DHOST}$ and $M_{500}^{\rm GR}$, plotted against the
redshift of the cluster. The horizontal dashed lines in the top and bottom
panels mark the GR case with $\Xi_{1}=0$ and $\tilde{\gamma}_{\rm N}=1$,
respectively. The shaded gray region represents a $10\%$ deviation from GR for
the mass parameter.

In Figure 2 we show the posterior probability density for the parameter
$\Xi_{1}$ obtained for each of the clusters. It is evident that the four non-
NFW clusters have posteriors notably far from the $\Xi_{1}=0$, with a $>95\%$
C.L upper limit of $\Xi_{1}<-0.5$. On the contrary, the other 8 clusters are
scattered around $\Xi_{1}=0$ showing no immediate preference for either a
negative or positive values for this parameter. Note that these results are
obtained under the prior assumption that NFW mass profile provides the best
description of data in the GR scenario, as shown in E19. Assuming an incorrect
or less-preferred mass profile lowers the quality of fit within the GR case
and could be partially compensated by the additional DHOST degree of freedom
showing a mislead preference for modifications to gravity. In this context, we
find it convenient that the validation for an assumption of the NFW mass
profile has already been performed in E19. It is very much possible that
modeling the four non-NFW clusters with a different profile would push the
constraints on $\Xi_{1}$ back towards the GR value. We do not perform this
analysis here since the large majority of the clusters we consider is well
fitted by the NFW profile. At a face-value the significance at which the non-
NFW mass profiles are preferred for these 4 clusters is much less than the
evidence for a DHOST modification we obtain here through the Bayesian evidence
(see Fig 2. of E19). However, we prefer to be conservative, and do not draw
any definitive conclusion from these considerations.

As described earlier, within the DHOST scenario, $M_{500}$ and
$\tilde{\gamma}_{\rm N}$ are degenerate parameters so that they are both not
well constrained when letting them free to independently vary in the the MCMC
analysis because of their expected strong correlation. Therefore, we also
present in Table 1 the derived constraints on $\tilde{\gamma}_{\rm N}\times
M_{500}$, which encompasses the total deviation from the GR case and will be
equal to the mass estimates in DHOST when $\tilde{\gamma}_{\rm N}=1$ is
assumed. As expected, one can immediately notice that $\tilde{\gamma}_{\rm
N}\times M_{500}$ is better constrained than the $M_{500}$ alone. In the
bottom panel of Figure 3, we show the redshift variation of the same,
normalized121212The normalization is performed by separately drawing $\sim
10000$ random samples from the MCMC chains of DHOST and the GR analyses.
Therefore, this is indicative of the total variation in $\tilde{\gamma}_{\rm
N}\times M_{500}$ from GR when marginalizing on $\Xi_{1}$, which is $0$ within
GR. to $M_{500}$ in the GR case. As in the case of $\Xi_{1}$, the four non-NFW
clusters once again show a large deviation from $\tilde{\gamma}_{\rm N}\times
M_{500}^{\rm DHOST}/M_{500}^{\rm GR}\sim 1$, which indicates either a very low
total mass in comparison to GR or that $\tilde{\gamma}_{\rm N}<<1$131313In
principle, the value of ${\gamma}_{\rm N}$ is expected to be $\sim 1$, unless
$1-\alpha_{\rm H}-3\beta_{1}<<1$, not implying a drastic change in $G_{\rm
N}^{\rm eff}$, which will be an outlier behavior. . The remaining 8 clusters,
do not immediately show any discernible trend in redsfhit for the estimates of
both $\Xi_{1}$ and $\tilde{\gamma}_{\rm N}\times M_{500}^{\rm
DHOST}/M_{500}^{\rm GR}$, however, see Section V.2 for an elaborate
discussion. A particular trend of these parameters in redshift could provide
strong implications for the redshift evolution of the physical parameters in
the DHOST scenario.

To assess the statistical evidence for/against the DHOST models, we compute
the Bayesian evidence ($\Delta\log{\cal B}$) with respect to the GR scenario.
These are reported in the last column of Table 1 clearly showing that GR is
mostly the preferred theory of gravity although DHOST is not strongly
disfavored. As for the constraints on $\Xi_{1}$, we still find that the
situation is reversed for the four non-NFW clusters. Note that the Bayesian
evidence is affected by the prior volume assumed on the parameters
$\Xi_{1},\,\tilde{\gamma}_{\rm N},\,M_{500}$ for the DHOST model. We,
therefore, perform a second MCMC analysis by fixing $\tilde{\gamma}_{\rm N}=1$
and sampling only over $\Xi_{1},\,M_{500}$ as free parameters. Firstly, we
notice that the posteriors for $\Xi_{1}$ remain equivalent to the ones
obtained with $\tilde{\gamma}_{\rm N}\neq 1$ for all the 12 clusters (see
e.g., Figure 5), indicating no correlation between the two effective
parameters of the DHOST model. We show the Bayesian evidence for the DHOST
analyses with $\tilde{\gamma}_{\rm N}\neq 1$ and the latter with
$\tilde{\gamma}_{\rm N}=1$ in Figure 4, indicated with filled and open
markers, respectively. While we see a small change/preference for the
$\tilde{\gamma}_{\rm N}=1$ w.r.t to the DHOST analysis with
$\tilde{\gamma}_{\rm N}\neq 1$, GR yet remains the preferred model of gravity.
Having established that there is no statistical preference for the DHOST
modification, from the individual X-COP clusters utilized in this work, we
proceed to discuss the nature of constraints for the DHOST modifications and
later a possible redshift dependence.

Figure 4: We show the Bayesian evidence estimates for the DHOST model w.r.t
GR. Here, DHOST analyses with $\tilde{\gamma}_{\rm N}\neq 1$ and
$\tilde{\gamma}_{\rm N}=1$ are indicated with filled and open markers,
respectively. To be contrasted using the Jeffery’s scale.

In Figure 5, we show the confidence regions for two clusters A2142 and
RXC1825, for the concentration and the DHOST parameters. Clearly, we have very
distinct distributions for the two clusters with RXC1825 showing more
degenerate boomerang-like contours in the parameter space of
$c_{500}\,\rm{vs.}\,\tilde{\gamma}_{\rm N}\times M_{500}$. We specifically
choose these two clusters A2142 and RXC1825, to represent $\Xi_{1}<0$ and
$\Xi_{1}>0$, cases (see also Appendix C), respectively. As can be clearly
noticed in Figure 5, the negative (positive) values of $\Xi_{1}$ indicate
strengthening (weakening) of gravity, which thereby implies a lower (higher)
value of $M_{500}^{\tilde{\gamma}_{\rm N}=1}$, with respect to the GR case. In
this context, we show the contours for the analysis with $\tilde{\gamma}_{\rm
N}=1$, essentially to validate that the constrains on the $\Xi_{1}$ parameter
remain unaltered, in comparison to the $\tilde{\gamma}_{\rm N}\neq 1$ case.
For comparison, we also show the constraints for the DHOST modification
obtained using only the $P_{\rm SZ}$ data. Through which we find that cluster
A2142 has more Gaussian-like constrains for the $M_{500}$, $\Xi_{1}$
parameters and hence is a simpler case for the discussion of extended gravity
theories141414Note that in (Ettori  _et al._ , 2017), a study of an emergent
gravity model was conducted for A2142 and A2319 clusters, as they are massive
clusters and makes them more suitable to study the modifications to gravity.
Which, however, was later extended to all the 12 clusters in X-COP in E19.,
while there is a shift in the constraints only for the $c_{500}$ parameter.
Note also that the constraints on $\tilde{\gamma}_{\rm N}\times M_{500}$
remain more consistent for all contrasted combinations, for both the clusters.
This we notice to be the case for all the clusters except the 4 non-NFW
clusters.

Using the $P_{\rm SZ}$ data alone, the RXC1825 cluster shows a multi-modal
behavior for the constrains on $c_{500}$ and consequently mild, yet a similar
double peak posterior for the $\Xi_{1}$ parameter. The posterior for the mass
however remains very much in agreement with the $P_{\rm SZ}+P_{\rm X}$
analysis, also for the GR case. This is also equivalently shown in Figure 9 of
Appendix D, where we show the distributions for $\Xi_{1}$ obtained for all the
clusters using $P_{\rm SZ}$ data alone. While posterior distribution on the
$\Xi_{1}$ appears as a double peak, it should be noted that the two peaks come
from very distinct distributions for the concentration parameter. The
inclusion of $P_{\rm X}$ data reduces the posteriors to the aforementioned
boomerang-like distributions, which highlights a change in the correlation of
the mass parameter with the $c_{500}$ and $\Xi_{1}$ parameters. Also
emphasizing the advantage of performing a joint analysis of $P_{\rm X}+P_{\rm
SZ}$ in our analysis. Interestingly, a similar effect was observed for the
degeneracy between the NFW parameters and the DHOST parameter $\Xi_{1}$ in a
simulated data based analysis in (Pizzuti  _et al._ , 2020)151515The analysis
therein is equivalent to the case of $\tilde{\gamma}_{\rm N}=1$, performed
here., however for much larger values of $\Xi_{1}$. In fact, this is the
reason for extended posteriors for $\Xi_{1}$ towards the positive end of the
distributions, as shown in Figure 2.

Also, it is interesting to note that the degeneracy manifests primarily
between the $c_{500}$ parameter and the DHOST parameter, rather than with the
$M_{500}$ parameter which remains consistent among different data combinations
within the DHOST analysis. In other words, $M_{500}$ and hence $R_{500}$ is
constrained in a more consistent way rather than the shape of the mass
profile, as it might be expected for the change in gravitational potential and
as an implication of a screening mechanism at play.

Figure 5:  We show the $68\%$ and $95\%$ C.L. contours for the concentration,
mass and DHOST parameters for the clusters A2142 (Left) and RXC1825 (Right).
Here $\cal{D}=\rm{I}$ stands for the dataset ($\cal{D}$) utilized being
$P_{\rm X}+P_{\rm SZ}$ and $\cal{D}=\rm{II}$, shows the constraints obtained
using only the $P_{\rm SZ}$ data. The dashed line marks the GR scenario of
$\Xi_{1}=0$.  Figure 6: Taylor expansion fit to the posterior distributions of
the $\Xi_{1}$ parameter. We show the 68$\%$ credible interval for the
$\Xi_{1}^{a}\neq 0$ (red) and $\Xi_{1}^{a}=0$ (blue) cases.

###  V.2 Redshift evolution of $\Xi_{1}$

While the analysis so far is done to constrain the DHOST parameters
independently for each of the clusters, they are clearly redshift dependent.
Alongside assessing the improvement of the fit to the pressure profile data
when including the DHOST modifications, it is also important to assess their
redshift behavior. We consider this fit to the posterior of $\Xi_{1}$ as the
first-order proxy of a possible redshift dependent behavior, which however
should be modeled simultaneously and is a computationally tedious analysis. In
this context, having a larger dataset and an extended range of redshift would
be of utmost importance to assess the DHOST gravity. Therefore, we fit the
posterior distributions on the $\Xi_{1}$ parameters, using the 8 clusters,
i.e., excluding the 4 non-NFW clusters to assess the redshift evolution of
$\Xi_{1}$. We consider a simple Taylor expansion as,

$\Xi_{1}(z)=\Xi_{1}^{0}+\frac{z}{1+z}\Xi_{1}^{\rm a},$ (21)

where $\Xi_{1}^{0}$ is value of $\Xi_{1}(z)$ at redshift $z=0$ and
$\Xi_{1}^{\rm a}\neq 0$ provides the redshift dependence. As can be seen from
Figure 6, there exists a redshift dependence, with $\Xi_{1}(z)$ deviating from
GR expectation for larger redshifts ($z\gtrsim 0.08$). However, note that the
redsfhit dependence is mostly dominated by cluster A2142, excluding which the
assumed functional form is consistent with no evolution within the redshift
range. In a similar approach fitting constant values through the 8 data
points, assuming no redsfhit dependence provides a joint
constrain161616Performing a joint analysis of all the 8 clusters with a
Gaussian prior on the $\Xi_{1}$ parameter would yield the same result.  on the
$\Xi_{1}$ parameter from the combination of the 8 clusters. We find this value
to be $\Xi_{1}=-0.030\pm 0.043$. This constraint is shown as the blue shaded
region in the Figure 6. Note that the fit presented here is an equivalent
approach as performing an importance sampling on the $\Xi_{1}$ parameter alone
constrained from the individual clusters. Which is a valid approach owing to
the fact that $\Xi_{1}$ is the only shared parameter among the clusters.

##  VI Comparison with existing constraints

It is worth noticing that the parameter $\Xi_{1}$ numerically coincides with
the $\Upsilon_{1}/4$ parameter in GLPV theories Sakstein  _et al._ (2016). The
form of the gravitational force related to the potential $\Phi$ in Equation 1
is unchanged, in the GLPV theory while the term with $\Xi_{3}$ is an
additional one. Therefore, for the observations where only the metric
potential $\Phi$ is involved, one can compare $\Xi_{1}$ with the existing
constraints in previous studies devoted to GLPV theory. For example, Koyama
and Sakstein (2015); Saito  _et al._ (2015) study the effects of modified
density profile for non-relativistic stars. On the other hand, astrophysical
observables, such as weak lensing that also depend on the potential $\Psi$
(see e.g. Koyama and Sakstein (2015); Sakstein  _et al._ (2016)), should be
taken into account in order to constrain $\Xi_{2}$ and $\Xi_{3}$ parameters.

Several works have earlier constrained the DHOST parameter $\Xi_{1}$, albeit
with varied parametric form171717For example, our parameter $\Xi_{1}$ here is
equivalent to $\Upsilon_{1}/4$ used in (Sakstein  _et al._ , 2016)., which
however we can readily compare with our results. As mentioned earlier in
Section I, in the current work we assess the modification to the hydrostatic
equilibrium alone and do not consider the weak lensing counterpart, that can
be utilized to place constraints on the $\Xi_{2}$ parameter. In this context,
(Sakstein  _et al._ , 2016), is one of the earliest analyses to provide
constraints on the $\Xi_{1},\,\Xi_{2}$ parameters by contrasting the
hydrostatic and weak lensing masses using a stacked dataset of 58 clusters in
the redshift range of $0.1\leq z\leq 1.2$. While the analysis in (Sakstein
_et al._ , 2016), improves the constraints by utilizing the stacked cluster
profiles, losing the redshift dependence, here we have for the first time (as
far as we are aware), tried to assess from individual unstacked clusters,
which helps us also assess the redshift dependence of the constraints. Our
constraints from individual clusters are however less stringent, and we find
good consistency with $\Xi_{1}\sim-0.028^{+0.23}_{-0.17}$, reported there
in181818Converted from the $\Upsilon_{1}=-0.11^{+0.98}_{-0.67}$ quoted in
(Sakstein  _et al._ , 2016).. The joint constraint of $\Xi_{1}=-0.030\pm
0.043$, derived from our proxy-fit to the posteriors of $\Xi_{1}$ from
individual clusters is in fact much more stringent, being $\sim 5$ times more
precise. This in turn, reasserts the promise in utilizing the clusters with a
combination of the X-ray data the SZ pressure profiles, which can improve upon
the existing constraints by almost an order in precision. Note also that the
redshift range ($z<0.1$) of the X-COP clusters used in our analysis is
complementary to the redshift range utilized in (Sakstein  _et al._ , 2016),
which is an added advantage when contrasting the constraints. Owing to the
mild differences in the mass estimates obtained in our analysis to those
quoted in the original analysis in E19, we also estimate a more conservative
limit by utilizing only 4 clusters (A1795, A3158, RXC1825 and ZW1215), which
agree very well. We find this limit to be $\Xi_{1}=-0.061\pm 0.074$, yet being
at least twice tighter than the previous estimate. It is worth noting that
these 4 clusters span a very small redshift range of $\Delta z\lesssim 0.02$
within $0.0597<z<0.0766$. For instance, the A1795 cluster provides an
individual constrain of $\Xi_{1}=-0.17^{+0.11}_{-0.10}$ using $P_{\rm
X}+P_{\rm SZ}$ data, which indicates a mild $\sim 1.55\,\sigma$ deviation from
GR, however statistically disfavored in comparison with GR with a Bayesian
evidence $\Delta\log(B)=-2.9$. As expected, we do not find such an evidence
constraint when using only the $P_{\rm X}$ data, as data from both the inner
and outer regions of the cluster need to be weighed against the model to
assess any variation in the gravitational potential along the radial
direction, while accounting for the screening effects.

Apart the cluster-scale constraints there exist lower limits from the non-
relativistic stars of $\Xi_{1}>-1/6$, requiring a stable static solution
(Saito  _et al._ , 2015) (see also (Saltas  _et al._ , 2018)). In fact, one
could reanalyze the current dataset by imposing the lower limit, which we
expect to, however, be consistent, except for the 4 non-NFW clusters. And an
upper limit of $\Xi_{1}<7\times 10^{-3}$ from the consistency of the minimum
mass for hydrogen burning in stars with the lowest mass red dwarf Sakstein
(2015a, b). Our constraints on the individual clusters (8 NFW) are very well
in accordance with these limits obtained from the much smaller scale objects.
Clusters A3266, RXC1825 and A85 have a positive mean value for $\Xi_{1}$ and
extend beyond the upper limit, however, being consistent with the values lower
than the limit as well. Taking the 4 non-NFW clusters at face value, one could
argue for a strong disagreement with the existing lower limit from the non-
relativistic stars, owing to i) non-NFW density profiles are not preferred as
strongly as DHOST modification from the Bayesian evidence (see Fig. 2. of E19)
and ii) they are not significantly affected by the non-thermal effects. This
would indicate invalidity of the aforementioned limits and scale-dependent
behavior of the current screening mechanism. However, we intend to interpret
these results more carefully for the moment, and wait to perform additional
analysis.

More recently, (Pizzuti  _et al._ , 2020) have performed a forecast analysis,
showing the potential galaxy cluster constrains for beyond Horndeski models,
scaling with the number of clusters. It has been shown that the posterior (
more precisely $\chi^{2}$ distribution ), could show large degeneracy for
large values of $\Xi_{1}$ and the NFW parameters used to model the mass
profile and that this should be taken into account when analyzing real
observations. However, in our analysis, we find the $\Xi_{1}$ parameter not
exceeding unity for any of the given 12 clusters. However, some confirming
correlations are seen in one or two clusters especially with the $P_{\rm SZ}$
data alone.

##  VII Summary

We constrain deviation from GR on the cluster scales, modeled through a
Vainshtein screening mechanism, utilizing the X-COP compilation which consists
of 12 galaxy clusters. While we have performed our analysis on all the 12
clusters, we make a careful selection of subsets consisting of either 8 or 4
clusters to comment on the final inferences: i) we eliminate 4 clusters as
they are not well represented by the NFW mass profile, ii) To be more
conservative, we exclude 4 other clusters which have mild $\sim 1\sigma$
variation from the masses quoted in the original analysis.

Our main results are summarized as follows:

  * •

Performing a Bayesian analysis using the so-called backward method assuming
the NFW mass profile, we find mild to moderate deviations from the GR scenario
($\Xi_{1}\neq 0$) of the order $\sim 2\sigma$ (see Table 1).

  * •

Comparing the Bayesian evidence, both GR and DHOST scenarios perform
equivalently, with mild to a moderate preference for GR at times.

  * •

As our main result, we report a constraint of $\Xi_{1}=-0.030\pm 0.043$
obtained using 8 clusters and a more conservative constraint of
$\Xi_{1}=-0.061\pm 0.074$ using only 4 clusters, for reasons quoted in the
text. Our result shows no indication for a deviation from GR, while being
stringent than the earlier constraints.

  * •

Assessing the redshift evolution we present the utility of datasets like the
one here. In effect, we find possible evidence for redshift dependent
behavior, however only dominated by constraints from A2142 cluster.

While individual clusters do not immediately suggest a significant
modification to GR, we find that a tentative redshift dependent behavior could
be observed at a larger significance. This indeed makes it essential that we
test the current setup against larger well-observed samples such as the NIKA2
SZ large program consisting of a sample size of 45 clusters, in a wider
redshift range. We also intend to extend the discussion in terms of the
physical parameters as outlined in Appendix B, which is of utmost importance
to assess the feasibility of the models within the current formalism.

## Acknowledgments

Authors are grateful to Stefano Ettori and Dominique Eckert for useful
comments and feedback on the use of data. BSH is supported by the INFN INDARK
grant. PK, MDP, VC, RM acknowledge support from INFN/Euclid Sezione di Roma.
PK, MDP and RM also acknowledge support from Sapienza Universitá di Roma
thanks to Progetti di Ricerca Medi 2018, prot. RM118164365E40D9 and 2019,
prot. RM11916B7540DD8D. We acknowledge the use of publicly available python
packages: numpy, scipy, & cmath.

## References

  * Adam  _et al._ (2016) R. Adam _et al._ (Planck), Astron. Astrophys. 594, A1 (2016), arXiv:1502.01582 [astro-ph.CO] .
  * Ade  _et al._ (2016a) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016a), arXiv:1502.01589 [astro-ph.CO] .
  * Riess  _et al._ (1998) A. G. Riess _et al._ (Supernova Search Team), Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ (1999) S. Perlmutter _et al._ (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * Bull  _et al._ (2016) P. Bull _et al._ , Phys. Dark Univ. 12, 56 (2016), arXiv:1512.05356 [astro-ph.CO] .
  * Caldwell  _et al._ (1998) R. R. Caldwell, R. Dave, and P. J. Steinhardt, Phys. Rev. Lett. 80, 1582 (1998), arXiv:astro-ph/9708069 .
  * Chiba  _et al._ (2000) T. Chiba, T. Okabe, and M. Yamaguchi, Phys. Rev. D 62, 023511 (2000), arXiv:astro-ph/9912463 .
  * Armendariz-Picon  _et al._ (2001) C. Armendariz-Picon, V. F. Mukhanov, and P. J. Steinhardt, Phys. Rev. D 63, 103510 (2001), arXiv:astro-ph/0006373 .
  * Tsujikawa (2010) S. Tsujikawa, Lect. Notes Phys. 800, 99 (2010), arXiv:1101.0191 [gr-qc] .
  * Nojiri and Odintsov (2011) S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59 (2011), arXiv:1011.0544 [gr-qc] .
  * Clifton  _et al._ (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Joyce  _et al._ (2016) A. Joyce, L. Lombriser, and F. Schmidt, Ann. Rev. Nucl. Part. Sci. 66, 95 (2016), arXiv:1601.06133 [astro-ph.CO] .
  * Langlois and Noui (2016) D. Langlois and K. Noui, JCAP 1602, 034 (2016), arXiv:1510.06930 [gr-qc] .
  * Crisostomi  _et al._ (2016) M. Crisostomi, K. Koyama, and G. Tasinato, JCAP 1604, 044 (2016), arXiv:1602.03119 [hep-th] .
  * Ben Achour  _et al._ (2016a) J. Ben Achour, D. Langlois, and K. Noui, Phys. Rev. D93, 124005 (2016a), arXiv:1602.08398 [gr-qc] .
  * Motohashi  _et al._ (2016) H. Motohashi, K. Noui, T. Suyama, M. Yamaguchi, and D. Langlois, JCAP 1607, 033 (2016), arXiv:1603.09355 [hep-th] .
  * Ben Achour  _et al._ (2016b) J. Ben Achour, M. Crisostomi, K. Koyama, D. Langlois, K. Noui, and G. Tasinato, JHEP 12, 100 (2016b), arXiv:1608.08135 [hep-th] .
  * Kobayashi (2019) T. Kobayashi, Rept. Prog. Phys. 82, 086901 (2019), arXiv:1901.07183 [gr-qc] .
  * Brans and Dicke (1961) C. Brans and R. H. Dicke, Phys. Rev. 124, 925 (1961).
  * De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Deffayet  _et al._ (2009) C. Deffayet, G. Esposito-Farese, and A. Vikman, Phys. Rev. D 79, 084003 (2009), arXiv:0901.1314 [hep-th] .
  * Neveu  _et al._ (2013) J. Neveu, V. Ruhlmann-Kleider, A. Conley, N. Palanque-Delabrouille, P. Astier, J. Guy, and E. Babichev, Astron. Astrophys. 555, A53 (2013), arXiv:1302.2786 [gr-qc] .
  * Neveu  _et al._ (2017) J. Neveu, V. Ruhlmann-Kleider, P. Astier, M. Besançon, J. Guy, A. Möller, and E. Babichev, Astron. Astrophys. 600, A40 (2017), arXiv:1605.02627 [gr-qc] .
  * Horndeski (1974) G. W. Horndeski, Int. J. Theor. Phys. 10, 363 (1974).
  * Kobayashi  _et al._ (2011) T. Kobayashi, M. Yamaguchi, and J. Yokoyama, Prog. Theor. Phys. 126, 511 (2011), arXiv:1105.5723 [hep-th] .
  * Zumalacárregui and García-Bellido (2014) M. Zumalacárregui and J. García-Bellido, Phys. Rev. D 89, 064046 (2014), arXiv:1308.4685 [gr-qc] .
  * Gleyzes  _et al._ (2015) J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, Phys. Rev. Lett. 114, 211101 (2015), arXiv:1404.6495 [hep-th] .
  * Vainshtein (1972) A. I. Vainshtein, Phys. Lett. B 39, 393 (1972).
  * Babichev and Deffayet (2013) E. Babichev and C. Deffayet, Class. Quant. Grav. 30, 184001 (2013), arXiv:1304.7240 [gr-qc] .
  * Brax  _et al._ (2004) P. Brax, C. van de Bruck, A.-C. Davis, J. Khoury, and A. Weltman, Phys. Rev. D 70, 123518 (2004), arXiv:astro-ph/0408415 .
  * Crisostomi and Koyama (2018) M. Crisostomi and K. Koyama, Phys. Rev. D97, 021301 (2018), arXiv:1711.06661 [astro-ph.CO] .
  * Langlois  _et al._ (2018a) D. Langlois, R. Saito, D. Yamauchi, and K. Noui, Phys. Rev. D 97, 061501 (2018a), arXiv:1711.07403 [gr-qc] .
  * Bartolo  _et al._ (2018) N. Bartolo, P. Karmakar, S. Matarrese, and M. Scomparin, JCAP 1805, 048 (2018), arXiv:1712.04002 [gr-qc] .
  * Dima and Vernizzi (2018) A. Dima and F. Vernizzi, Phys. Rev. D97, 101302 (2018), arXiv:1712.04731 [gr-qc] .
  * Hirano  _et al._ (2019) S. Hirano, T. Kobayashi, and D. Yamauchi, Phys. Rev. D 99, 104073 (2019), arXiv:1903.08399 [gr-qc] .
  * Crisostomi  _et al._ (2019) M. Crisostomi, M. Lewandowski, and F. Vernizzi, Phys. Rev. D 100, 024025 (2019), arXiv:1903.11591 [gr-qc] .
  * Kimura  _et al._ (2012) R. Kimura, T. Kobayashi, and K. Yamamoto, Phys. Rev. D 85, 024023 (2012), arXiv:1111.6749 [astro-ph.CO] .
  * Narikawa  _et al._ (2013) T. Narikawa, T. Kobayashi, D. Yamauchi, and R. Saito, Phys. Rev. D 87, 124006 (2013), arXiv:1302.2311 [astro-ph.CO] .
  * Koyama  _et al._ (2013) K. Koyama, G. Niz, and G. Tasinato, Phys. Rev. D 88, 021502 (2013), arXiv:1305.0279 [hep-th] .
  * Kase and Tsujikawa (2013) R. Kase and S. Tsujikawa, JCAP 08, 054 (2013), arXiv:1306.6401 [gr-qc] .
  * Kobayashi  _et al._ (2015) T. Kobayashi, Y. Watanabe, and D. Yamauchi, Phys. Rev. D 91, 064013 (2015), arXiv:1411.4130 [gr-qc] .
  * Koyama and Sakstein (2015) K. Koyama and J. Sakstein, Phys. Rev. D91, 124066 (2015), arXiv:1502.06872 [astro-ph.CO] .
  * Saito  _et al._ (2015) R. Saito, D. Yamauchi, S. Mizuno, J. Gleyzes, and D. Langlois, JCAP 1506, 008 (2015), arXiv:1503.01448 [gr-qc] .
  * Sakstein (2015a) J. Sakstein, Phys. Rev. Lett. 115, 201101 (2015a), arXiv:1510.05964 [astro-ph.CO] .
  * Sakstein (2015b) J. Sakstein, Phys. Rev. D92, 124045 (2015b), arXiv:1511.01685 [astro-ph.CO] .
  * Jain  _et al._ (2016) R. K. Jain, C. Kouvaris, and N. G. Nielsen, Phys. Rev. Lett. 116, 151103 (2016), arXiv:1512.05946 [astro-ph.CO] .
  * Babichev  _et al._ (2016) E. Babichev, K. Koyama, D. Langlois, R. Saito, and J. Sakstein, Class. Quant. Grav. 33, 235014 (2016), arXiv:1606.06627 [gr-qc] .
  * Sakstein  _et al._ (2016) J. Sakstein, H. Wilcox, D. Bacon, K. Koyama, and R. C. Nichol, JCAP 1607, 019 (2016), arXiv:1603.06368 [astro-ph.CO] .
  * Salzano  _et al._ (2017) V. Salzano, D. F. Mota, S. Capozziello, and M. Donahue, Phys. Rev. D 95, 044038 (2017), arXiv:1701.03517 [astro-ph.CO] .
  * Eckert  _et al._ (2017) D. Eckert, S. Ettori, E. Pointecouteau, S. Molendi, S. Paltani, and C. Tchernin, Astron. Nachr. 338, 293 (2017), arXiv:1611.05051 [astro-ph.CO] .
  * Ettori  _et al._ (2019) S. Ettori, V. Ghirardini, D. Eckert, E. Pointecouteau, F. Gastaldello, M. Sereno, M. Gaspari, S. Ghizzardi, M. Roncarelli, and M. Rossetti, Astron. Astrophys. 621, A39 (2019), arXiv:1805.00035 [astro-ph.CO] .
  * Eckert  _et al._ (2019) D. Eckert _et al._ , Astron. Astrophys. 621, A40 (2019), arXiv:1805.00034 [astro-ph.CO] .
  * Ghirardini  _et al._ (2019) V. Ghirardini _et al._ , Astron. Astrophys. 621, A41 (2019), arXiv:1805.00042 [astro-ph.CO] .
  * Ade  _et al._ (2016b) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A24 (2016b), arXiv:1502.01597 [astro-ph.CO] .
  * Terukina  _et al._ (2014) A. Terukina, L. Lombriser, K. Yamamoto, D. Bacon, K. Koyama, and R. C. Nichol, JCAP 04, 013 (2014), arXiv:1312.5083 [astro-ph.CO] .
  * Wilcox  _et al._ (2015) H. Wilcox _et al._ , Mon. Not. Roy. Astron. Soc. 452, 1171 (2015), arXiv:1504.03937 [astro-ph.CO] .
  * Pizzuti  _et al._ (2020) L. Pizzuti, I. D. Saltas, and L. Amendola, (2020), arXiv:2011.15089 [astro-ph.CO] .
  * Langlois  _et al._ (2018b) D. Langlois, R. Saito, D. Yamauchi, and K. Noui, Phys. Rev. D97, 061501 (2018b), arXiv:1711.07403 [gr-qc] .
  * Langlois  _et al._ (2017) D. Langlois, M. Mancarella, K. Noui, and F. Vernizzi, JCAP 1705, 033 (2017), arXiv:1703.03797 [hep-th] .
  * Kravtsov and Borgani (2012) A. V. Kravtsov and S. Borgani, ARA&A 50, 353 (2012), arXiv:1205.5556 [astro-ph.CO] .
  * Vikhlinin  _et al._ (2006) A. Vikhlinin, A. Kravtsov, W. Forman, C. Jones, M. Markevitch, S. S. Murray, and L. Van Speybroeck, Astrophys. J. 640, 691 (2006), arXiv:astro-ph/0507092 [astro-ph] .
  * Vikhlinin (2006) A. Vikhlinin, Astrophys. J. 640, 710 (2006), arXiv:astro-ph/0504098 [astro-ph] .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Ade  _et al._ (2014) P. Ade _et al._ (Planck), Astron. Astrophys. 571, A29 (2014), arXiv:1303.5089 [astro-ph.CO] .
  * Ghirardini  _et al._ (2018) V. Ghirardini, S. Ettori, D. Eckert, S. Molendi, F. Gastaldello, E. Pointecouteau, G. Hurier, and H. Bourdin, Astron. Astrophys. 614, A7 (2018), arXiv:1708.02954 [astro-ph.CO] .
  * Eckert  _et al._ (2016) D. Eckert _et al._ , Astron. Astrophys. 592, A12 (2016), arXiv:1512.03814 [astro-ph.CO] .
  * Eckert  _et al._ (2015) D. Eckert, S. Ettori, S. Molendi, F. Vazza, M. Roncarelli, F. Gastaldello, and M. Rossetti, Mon. Not. Roy. Astron. Soc. 447, 2198 (2015), arXiv:1310.8389 [astro-ph.CO] .
  * Hogg and Foreman-Mackey (2018) D. W. Hogg and D. Foreman-Mackey, Astrophys. J. Suppl. 236, 11 (2018), arXiv:1710.06068 [astro-ph.IM] .
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Hinton (2016) S. R. Hinton, The Journal of Open Source Software 1, 00045 (2016).
  * Trotta (2008) R. Trotta, Contemp. Phys. 49, 71 (2008), arXiv:0803.4089 [astro-ph] .
  * Trotta (2017) R. Trotta (2017) arXiv:1701.01467 [astro-ph.CO] .
  * Heavens  _et al._ (2017a) A. Heavens, Y. Fantaye, E. Sellentin, H. Eggers, Z. Hosenie, S. Kroon, and A. Mootoovaloo, Phys. Rev. Lett. 119, 101301 (2017a), arXiv:1704.03467 [astro-ph.CO] .
  * Heavens  _et al._ (2017b) A. Heavens, Y. Fantaye, A. Mootoovaloo, H. Eggers, Z. Hosenie, S. Kroon, and E. Sellentin, (2017b), arXiv:1704.03472 [stat.CO] .
  * Jeffreys (1961) H. Jeffreys,  _The Theory of Probability_ , Oxford Classic Texts in the Physical Sciences (1961).
  * Ameglio  _et al._ (2009) S. Ameglio, S. Borgani, E. Pierpaoli, K. Dolag, S. Ettori, and A. Morandi, Mon. Not. Roy. Astron. Soc. 394, 479 (2009), arXiv:0811.2199 [astro-ph] .
  * Croston  _et al._ (2006) J. H. Croston, M. Arnaud, E. Pointecouteau, and G. W. Pratt, Astron. Astrophys. 459, 1007 (2006), arXiv:astro-ph/0608700 .
  * Mathiesen  _et al._ (1999) B. Mathiesen, A. E. Evrard, and J. J. Mohr, Astrophys. J. Lett. 520, L21 (1999), arXiv:astro-ph/9904429 .
  * Ettori  _et al._ (2013) S. Ettori, A. Donnarumma, E. Pointecouteau, T. Reiprich, S. Giodini, L. Lovisari, and R. Schmidt, Space Sci. Rev. 177, 119 (2013), arXiv:1303.3530 [astro-ph.CO] .
  * Ettori  _et al._ (2017) S. Ettori, V. Ghirardini, D. Eckert, F. Dubath, and E. Pointecouteau, Mon. Not. Roy. Astron. Soc. 470, L29 (2017), arXiv:1612.07288 [astro-ph.CO] .
  * Saltas  _et al._ (2018) I. D. Saltas, I. Sawicki, and I. Lopes, JCAP 05, 028 (2018), arXiv:1803.00541 [astro-ph.CO] .
  * Arai  _et al._ (2020) S. Arai, P. Karmakar, and A. Nishizawa, Phys. Rev. D 102, 024003 (2020), arXiv:1912.01768 [gr-qc] .

##  Appendix A DHOST action

The action for the viable Class Ia DHOST theory after GW170817 event
($c_{g}^{2}=c^{2}$) is

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
d^{4}x\,\sqrt{-g}\,\mathcal{L}\,,$ (22)

where

$\begin{split}&L^{{}_{\rm
DHOST}}_{c_{g}^{2}=c^{2}}=P+Q\,\Box\phi+F\,R+A_{3}\phi^{\mu}\phi^{\nu}\phi_{\mu\nu}\Box\phi\\\
&+\frac{1}{8F}\bigg{(}48{F_{X}}^{2}-8(F-X{F_{X}})A_{3}-X^{2}A_{3}^{2}\bigg{)}\phi^{\mu}\phi_{\mu\nu}\phi_{\lambda}\phi^{\lambda\nu}\\\
&+\frac{1}{2F}\left(4{F_{X}}+XA_{3}\right)A_{3}(\phi_{\mu}\phi^{\mu\nu}\phi_{\nu})^{2}\;.\end{split}$
(23)

$P,Q,F,A_{3}$ are the arbitrary functions of the scalar field $\phi$ and its
kinetic energy, $X$, which reduces to the GR case when we set $F=1/2\kappa$
with $\kappa=8\pi G_{\rm N}c^{-4}$, and $P,Q,A_{3}=0$.

In effective field theory (EFT) formalism, the action is expressed in terms of
time dependent linear operators. $c_{g}=c$ tightly constrains the tensor speed
alteration parameter, $\alpha_{\rm{T}}$. The viable DHOST theory ($c_{g}=c$)
is expressed in terms of five time-dependent linear effective field theory
(EFT) parameters, i.e., $\alpha_{\rm{M,B,K,H}}$ and $\beta_{1}$, which are the
measures to the deviation from the $\Lambda$CDM model Langlois  _et al._
(2017). The coefficient of the linear perturbations only depends on the
background value, hence only on the time. Subsequently,
$\alpha_{\rm{M,B,K,H}}$ and $\beta_{1}$ are solely functions of time.

##  Appendix B A note on change of basis

Figure 7: We show the scatter for the physical parameters of the DHOST theory,
$\\{\alpha_{\rm H},\,\beta_{1}\\}$. See Appendix B for discussion.

While we have performed the analysis using the effective parameters
$\\{\Xi_{1},\tilde{\gamma}_{\rm N}\\}$, they are both functions of the
physical parameters $\\{\alpha_{\rm H},\,\beta_{1},\,{\gamma}_{\rm N}\\}$,
which in turn are functions of time, and hence redshift. As mentioned earlier,
$\tilde{\gamma}_{\rm N}$ is unconstrained when sampled over as a free
parameter and the posterior for $\Xi_{1}$ remains unchanged when fixing
$\tilde{\gamma}_{\rm N}=1$. Therefore one could advantage of the analysis
performed assuming $\tilde{\gamma}_{\rm N}=1$, and estimate,

$\tilde{\gamma}_{\rm N}\equiv\frac{\tilde{\gamma}_{\rm N}\times
M_{500}^{\tilde{\gamma}_{\rm N}\neq 1}}{M_{500}^{\tilde{\gamma}_{\rm N}=1}}.$
(24)

Doing so, we find strictly a mean value of $\tilde{\gamma}_{\rm N}\sim 1$ as
expected, however with the dispersion that varies from cluster to cluster. If
one were to assume a fixed but $\tilde{\gamma}_{\rm N}\neq 1$ value, the mean
obtained from the above expression would be expected to be same as the
assumption. Therefore, we infer this dispersion as the tentative uncertainty
on the GR expectation of $\tilde{\gamma}_{\rm N}=1$, in contrast to having
either a fixed value or a completely unconstrained quantity. This validates
the fact that the $\Xi_{1}$ and $\tilde{\gamma}_{\rm N}$ are uncorrelated
parameters and hence reasserts the method followed using two different MCMC
samples to estimate this dispersion. In other words, it is an equivalent
approach assuming that the mass of the cluster is accurately known in the
DHOST gravity and that the uncertainty is present only on the
$\tilde{\gamma}_{\rm N}$ parameter. Needless to say, having probes such as
weak lensing or any background expansion history, could break the degeneracy
and help correctly constrain this quantity. Here we only anticipate an
analysis assigning the uncertainty on the parameter $\tilde{\gamma}_{\rm
N}\times M_{500}$ completely to $\tilde{\gamma}_{\rm N}$ alone, however to
obtain an expectation for the allowed parameter space of the physical
quantities ($\alpha_{\rm H},\,\beta_{1}$) within the DHOST theory.

Through a simple change of basis we can represent the $\alpha_{\rm H}$ and
$\beta_{1}$ in terms of effective parameters $\\{\Xi_{1},\tilde{\gamma}_{\rm
N}\\}$ as,

$\displaystyle\begin{aligned} \alpha_{\rm
H}&=\frac{1}{2}\left[\left(-1+\frac{{\gamma}_{\rm N}}{\tilde{\gamma}_{\rm
N}}+\frac{3\Xi_{1}}{2}\right)\pm\frac{3}{2}\sqrt{\xi}\right]\\\
\beta_{1}&=\frac{1}{2}\left[\left(1-\frac{{\gamma}_{\rm
N}}{\tilde{\gamma}_{\rm
N}}-\frac{\Xi_{1}}{2}\right)\mp\frac{1}{2}\sqrt{\xi}\right]\\\
\textrm{where,}\\\ \xi&=\Xi_{1}^{2}+4\Xi_{1}\left(1-\frac{{\gamma}_{\rm
N}}{\tilde{\gamma}_{\rm N}}\right).\end{aligned}$

Note that the above solutions for $\alpha_{\rm H}$ and $\beta_{1}$ are valid
under the conditions, $\left[\Xi_{1}>0,\,\,\tilde{\gamma}_{\rm
N}>{\gamma}_{\rm N}\times 4/(4+\Xi_{1})\right]$ or
$\left[\Xi_{1}<0,\,\,\tilde{\gamma}_{\rm N}<{\gamma}_{\rm N}\times
4/(4+\Xi_{1})\right]$. While we have the distribution on $\tilde{\gamma}_{\rm
N}$ and $\Xi_{1}$, the parameter ${\gamma}_{\rm N}$ itself is free quantity
and therefore certain values for the same have to be assumed before obtaining
the distributions of $\\{\alpha_{\rm H},\beta_{1}\\}$. As in, the scatter for
values of $\\{\alpha_{\rm H},\,\beta_{1}\\}$, for viable physical models can
obtained using a flat distribution of $0.8<{\gamma}_{\rm N}<1.2$. Note that
the phenomenology of DHOST gravity can be equivalent to that of GR, even if
$\alpha_{\rm H},\,\beta_{1}\neq 0$, when $\alpha_{\rm H}=-\beta_{1}$ and
${\gamma}_{\rm N}=1-2\beta_{1}\equiv 1+2\alpha_{\rm H}$. The distributions of
DHOST parameters have been estimated by computing the approximated time
evolution of a scalar field through numerical simulations in Arai  _et al._
(2020). A comparison can be done with simulations and data from future
observations. We leave the possibility of a detailed analysis of the physical
parameters along these lines for a future consideration, also with the
inclusion of complementary weak lensing probes, which would allow to break
additional degeneracies.

Figure 8: Here we show the comparison of the mass estimates ($M_{500}$),
obtained in our analysis for both GR and DHOST cases, against the GR estimates
in E19.

##  Appendix C Mass comparison

In this section we provide a comparison the $M_{500}$ estimates obtained
within our analysis and the original analysis presented in E19, as in Figure
8. We also compare constrains in the DHOST scenario, which are shown as open
markers. Firstly we notice an overall agreement for all the clusters in the GR
case. However interestingly, we find that in our analysis, the uncertainties
estimated for the non-NFW clusters are lower, while the mean values are in
better agreement. Also the DHOST mass estimates for these 4 clusters tend to
follow a constant scaling w.r.t the GR masses, of $M_{500}^{\rm DHOST}\sim
0.6\times M_{500}^{\rm GR}$. The effects of mass variation are also
represented equivalently in the lower panel of Figure 3 as a ratio. This is
indeed an interesting feature, which requires more attention and we intend to
explore this in a future study.

The upper panels in Figure 10, show the mass reconstructions both in the GR
(blue) and the DHOST (red) scenarios. In here one can notice that mass
reconstructions in the DHOST scenario for the non-NFW clusters show larger
variation from the GR case, and this illustrates the large values of the
Bayesian evidence in favor of the DHOST modeling, reported in Table 1 (see
also Figure 4). Th lower panels, also illustrate how the constraining ability
of the $P_{\rm SZ}$ and $P_{\rm X}$ data independently contribute to the joint
likelihood (see Equation 19), through the variation of uncertainty in the
radial distance.

As mentioned earlier we have utilized the Vikhlinin profile with only 6
parameters neglecting a second component (see eq.3 of (Vikhlinin  _et al._ ,
2006)). Firstly, we find that the agreement between our mass estimates here
and the mass estimates reported in E19 worsens when we utilize the full
parametric form with 9 free parameters. Alongside this, we also obtain much
larger uncertainties on the mass estimates, which deters us from further using
the full parametric form in the main analysis.

Figure 9: Same as Figure 2, using only the $P_{\rm SZ}$ data. Note here that
the limits on the x-axis are different from those in Figure 2. Also, we have
not excluded the three inner most data points as done in the main analysis.

##  Appendix D Constraints using $P_{\rm SZ}$ data alone

Alongside our main analysis, we also assess the constraints on the DHOST
parameter $\Xi_{1}$ when utilizing the $P_{\rm SZ}$ data alone. In Figure 9,
we show the marginalized posterior distributions, which can be compared with
the main results in Figure 2. In here, however, we include also the 3 internal
points of the $P_{\rm SZ}$ data, which were excluded in the main analysis.
Notice that the constraints on the parameter $\Xi_{1}$, vary significantly and
we find good consistency between our main analysis and the $P_{\rm SZ}$
analysis only for the clusters, A2142, RXC1825. For the rest of the clusters,
however, the constraints vary, signifying the need for the formalism followed
in the main analysis, i.e., exclusion of 3 innermost points of $P_{\rm SZ}$
data and joint analysis of $P_{\rm SZ}$ and $P_{\rm X}$. We no longer notice
the significant clustering of the constraints on $\Xi_{1}\ll 0$ obtained for
the non-NFW clusters and in general the constraints from individual clusters
are more dispersed.

Figure 10: We show the comparison of mass profiles (Top) and the relative
difference (Bottom), between the GR and the DHOST modification. The dashed and
dash-dotted vertical lines show the $R_{500}$ in the GR and DHOST cases,
respectively. The shaded region corresponds to a $68\%$ C.L. limits on the
mass profiles. The corresponding mass ($M_{500}$) and concentration
($c_{500}$) parameters are shown in Table 1.