#  Caustic and hydrostatic mass bias: Implications for modified gravity

Minahil Adil Butt<EMAIL_ADDRESS>SISSA-International School for Advanced
Studies, Via Bonomea 265, 34136 Trieste, Italy  ICTP-The Abdus Salam
International Centre for Theoretical Physics, Strada Costiera 11, 34151
Trieste, Italy  IFPU, Institute for Fundamental Physics of the Universe, via
Beirut 2, 34151 Trieste, Italy  INFN, Sezione di Trieste, Via Valerio 2,
I-34127 Trieste, Italy  Balakrishna S. Haridasu<EMAIL_ADDRESS>SISSA-
International School for Advanced Studies, Via Bonomea 265, 34136 Trieste,
Italy  IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2,
34151 Trieste, Italy  INFN, Sezione di Trieste, Via Valerio 2, I-34127
Trieste, Italy  Yacer Boumechta  SISSA-International School for Advanced
Studies, Via Bonomea 265, 34136 Trieste, Italy  ICTP-The Abdus Salam
International Centre for Theoretical Physics, Strada Costiera 11, 34151
Trieste, Italy  IFPU, Institute for Fundamental Physics of the Universe, via
Beirut 2, 34151 Trieste, Italy  INFN, Sezione di Trieste, Via Valerio 2,
I-34127 Trieste, Italy  Francesco Benetti  SISSA-International School for
Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy  Lorenzo Pizzuti
<EMAIL_ADDRESS>Dipartimento di Fisica G. Occhialini, Universitá degli Studi
di Milano Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy  Carlo
Baccigalupi<EMAIL_ADDRESS>SISSA-International School for Advanced Studies,
Via Bonomea 265, 34136 Trieste, Italy  IFPU, Institute for Fundamental Physics
of the Universe, via Beirut 2, 34151 Trieste, Italy  INFN, Sezione di Trieste,
Via Valerio 2, I-34127 Trieste, Italy  Andrea Lapi<EMAIL_ADDRESS>SISSA-
International School for Advanced Studies, Via Bonomea 265, 34136 Trieste,
Italy  IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2,
34151 Trieste, Italy  INFN, Sezione di Trieste, Via Valerio 2, I-34127
Trieste, Italy  IRA-INAF, Via Gobetti 101, 40129 Bologna, Italy

###### Abstract

We propose and perform a joint analysis of the two different mass estimates of
galaxy clusters, namely the hydrostatic and caustic techniques. Firstly, we
show comprehensively that the mass bias between these two techniques can be
possibly alleviated when cluster-specific assumptions constrained using the
hydrostatic technique are utilized within the caustic technique. While at face
value this demotes the caustic technique from a completely independent method,
this allows one to further tighten the constraints on the cluster mass and
subsequently, allow us to test modifications to gravity. Implementing the
aforementioned formalism for two well-observed massive galaxy clusters, A2029
and A2142, we highlight the proof of concept. In the current implementation,
we use this method to constrain the Chameleon screening and Vainshtein
screening. As anticipated, we show that the joint analysis can help improve
the constraints on these modified gravity scenarios.

##  I Introduction

The accurate determination of galaxy cluster masses stands as a cornerstone in
our quest to understand the evolution of cosmic structures and exploit
clusters Eckert  _et al._ (2022); Aarseth  _et al._ (1979); de Andres  _et
al._ (2022) as indispensable cosmological tools Cui  _et al._ (2018); Ade  _et
al._ (2016); Adam  _et al._ (2016). Beyond their significance in the broader
cosmological context, cluster mass determination holds a pivotal role in
deciphering the intricate astrophysical processes governing both baryonic and
dark matter physics within these clusters, encompassing the hot gas of the
intracluster medium (ICM) and member galaxies, as well as the thermodynamical
properties of the clusters in the outskirts Eckert  _et al._ (2015, 2016);
Ettori  _et al._ (2013, 2019a).

Cluster mass measurement, however, presents an intricate challenge, primarily
because the dominant dark matter component can only be indirectly probed,
often relying on assumed fitting functions such as the Navarro-Frenk-White
(NFW) profile Navarro  _et al._ (1996). Determining the total mass of a
cluster traditionally hinges upon the gravitational influence it exerts on the
properties of the ICM and member galaxies Ettori  _et al._ (2019b), or its
gravitational lensing effect on the light originating from background sources
Murray  _et al._ (2022). The most accurate and precise mass estimation
techniques include hydrostatic masses determined from X-ray observations of
the ICM, caustic techniques based on galaxy dynamics Diaferio (1999),
kinematics of the galaxy cluster members solving the Jeans equation (for
example Biviano  _et al._ (2023)) and weak gravitational lensing measurements
von der Linden  _et al._ (2014).

The gold standard for cluster mass estimation has long been hydrostatic masses
derived from X-ray observations of the ICM Ettori  _et al._ (2019b).
Nevertheless, the accuracy of these estimates can be compromised by departures
from hydrostatic equilibrium or the presence of non-thermal pressure sources,
including turbulence and cosmic rays, often leading to systematic
underestimations of the true mass by 10-30 percent Ettori  _et al._ (2019b);
Hurier and Angulo (2018). Recent advances in both techniques and data quality
have spotlighted these biases, sparking a reevaluation of hydrostatic mass
estimates against alternative methods such as weak gravitational lensing
Maughan  _et al._ (2016); Ettori  _et al._ (2019b). Among the most robust mass
estimation techniques, weak gravitational lensing measurements, which probe
the distortion of background source light due to the cluster’s gravitational
field, offer unparalleled precision and independence from dynamical
assumptions. Nevertheless, they too necessitate high-quality data, including
extensive galaxy redshifts and lensed sources, as well as assumptions in de-
projection to reconstruct the shape of the three-dimensional mass profile
(e.g. Umetsu (2020)).
This paper embarks on an exploration of caustic techniques for mass profiling,
serving as a promising complement to weak gravitational lensing et al (2018);
et al. (2018). The caustic method relies on identifying unique structures
within line-of-sight velocity and projected-radius space that mirror the
cluster’s escape velocity profile. These caustics provide a means to
reconstruct the cluster’s enclosed mass, extending far beyond the virial
radius Diaferio (1999); Gifford  _et al._ (2013); Serra  _et al._ (2010).
Similar to the lensing measurements et al. (2023, 2012); Diaferio  _et al._
(2005), caustic masses remain agnostic to the dynamical state of clusters and
are resilient against the physical processes that can bias hydrostatic mass
estimates. In practice, the caustic technique capitalizes on overdense
envelopes in phase-space diagrams, mapping the line-of-sight velocity versus
projected radius for galaxies inhabiting cluster infall regions Diaferio and
Geller (1997); Diaferio (1999, 2009). These envelopes’ edges, referred to as
the caustic profile ($\mathcal{A}(r)$), can be extracted from observational
data Serra  _et al._ (2010), though needing intricate techniques. An inherent
assumption in caustic mass estimates is the constancy of the filling factor
$\mathcal{F}_{\beta}(r)$, a parameter linked to the ratio of mass gradient to
gravitational potential Diaferio (1999); Serra  _et al._ (2010). Nevertheless,
N-body simulations have unveiled breakdowns in this approximation within
cluster inner regions, potentially resulting in an overestimate of up to 10-20
percent at $r_{500}$, i.e. the radius enclosing an average overdensity of
$500$ times the critical density of the universe at that redshift (see Section
V), increasing at smaller radii Serra  _et al._ (2010). While hydrostatic and
caustic methods operate independently, each resting on distinct assumptions
and systematic uncertainties, they exhibit reasonable agreement on average
Logan  _et al._ (2022); Andreon  _et al._ (2017). While a significant scatter
is present between the two estimators, mass estimates align to within
approximately 20 percent over the full radial range sampled by both techniques
Maughan  _et al._ (2016).

In the pursuit of alleviating biases between hydrostatic and caustic masses
Maughan  _et al._ (2016); Logan  _et al._ (2022), we speculate on the effects
of relaxing the constant $\mathcal{F}_{\beta}(r)$ assumption Diaferio (1999);
Gifford  _et al._ (2013). To this end, we explore varying
$\mathcal{F}_{\beta}(r)$ under the framework of General Relativity (GR)
(assuming a $\Lambda$CDM Universe), and two modified gravity models: Chameleon
Terukina  _et al._ (2014) and Vainshtein Screening Babichev and Deffayet
(2013). The Lagrangian of the Chameleon screening theory includes the usual
Einstein-Hilbert Lagrangian plus the scalar field, in addition to the Standard
Model fields coupled minimally to gravity Khoury and Weltman (2004a);
Zaregonbadi  _et al._ (2022); Ivanov and Wellenzohn (2020); Kraiselburd  _et
al._ (2018); Tsujikawa  _et al._ (2009). In this model we obtain interesting
phenomenology on cosmological scales while simultaneously screening MG effects
in high-density environments, allowing it to avoid solar system constraints.
Thus, GR is obtained in regions of high-density while in low-density
environments we observe the effects of the modified potentials. The Vainshtein
screening Vainshtein (1972) is a useful mechanism for the higher-order scalar
tensor Babichev and Deffayet (2013); Brax  _et al._ (2004). In Vainshtein
screening, the gravitational potentials are modified inside the matter sources
Crisostomi and Koyama (2018); Langlois  _et al._ (2018); Bartolo  _et al._
(2018); Dima and Vernizzi (2018); Hirano  _et al._ (2019); Crisostomi  _et
al._ (2019). This mechanism allows to hide via non-linear effects - typically
for source distances smaller than a so-called Vainshtein radius which depends
on the source and the theory considered - some degrees of freedom whose
effects are then only left important at large distances, e.g. for cosmology.
Both these theories have been tested using galaxy clusters Terukina  _et al._
(2014); Haridasu  _et al._ (2023).

Our analysis aims to shed light on the resultant impact on the caustic mass
and, by extension, the mass bias, scrutinizing diverse regions of the modified
gravity parameter space. The structure of this paper unfolds as follows:
Section II offers a concise overview of the Caustic technique. In Section III,
we delve into the intricacies of screening mechanisms within the framework of
modified gravity scenarios. Section V outlines our approach, encompassing data
sources, methodology, and analytical techniques. Our results, from the current
investigation, are presented in Section VI. Finally, Section VII provides
concluding remarks. Throughout this work, we assume $H_{0}=70\,\text{km/s
Mpc${}^{-1}$}$ and $\Omega_{m}=0.3$ in the estimation of the critical density
$\rho_{\rm{crit}}=8\pi G/3H(z)^{2}$ with the standard $\Lambda$CDM background
having $H(z)^{2}=\Omega_{m}(1+z)^{3}+1-\Omega_{m}$.

##  II The Caustic Method

The caustic technique, as introduced by Diaferio and Geller (1997) and further
developed in subsequent works (Diaferio (1999, 2009); Serra  _et al._ (2010)),
provides a unique approach to estimate the mass of the clusters, utilizing the
escape velocity profile of member galaxies within the cluster. This technique
extends our insights from the central cluster region to radii as large as
three times the cluster’s virial radius, denoted as $r_{200}$, where $r_{200}$
signifies the radius of a sphere with an average density 200 times the
critical density ($\rho_{\rm{crit}}$) of the background universe. This is
particularly valuable because, at such large radii, galaxy clusters might not
be in complete dynamical equilibrium, based solely on galaxy redshift data.
Hierarchical clustering, a prominent formation mechanism for galaxy clusters,
involves the aggregation of smaller systems. Unlike a purely radial infall
expected in the spherical collapse model, this process incorporates
substantial non-radial velocities (e.g., White et al. 2010), resulting in a
complex velocity distribution of the galaxies within clusters (see reviews in
Diaferio (2009) and Serra  _et al._ (2010)). The primary reason for this
distribution is completely dependent on the local gravitational potential
Diaferio and Geller (1997), coupled with the influences of surrounding groups
and tidal fields. When visualized in a redshift diagram, portraying the line-
of-sight velocity versus projected distance from the cluster center (the so-
called projected phase space, p.p.s.), cluster members delineate a distinctive
trumpet-shaped region symmetric along the radial axis Kaiser (1987); van
Haarlem and van de Weygaert (1993); Regos and Geller (1989).
Within this framework, the caustic surface serves as a crucial demarcation,
defining the boundaries of this trumpet-shaped region. The amplitude of these
caustics, denoted as $\mathcal{A}(r)$, diminishes as we move away from the
cluster center, and it is intrinsically tied to the average velocity component
$\langle v^{2}\rangle$ Diaferio  _et al._ (2005). In a spherically symmetric
system, the escape velocity $v_{\rm{esc}}^{2}(r)$ is directly related to the
gravitational potential $\Phi(r)$ and is a non-increasing function of radial
distance. The highest observable velocity corresponds to the escape velocity.
Therefore, $\mathcal{A}^{2}(r)$ at a projected radius $r_{\perp}$ effectively
measures the average velocity component along the line of sight at the three-
dimensional radius $r=r_{\perp}$.

To quantify this average velocity component, the velocity anisotropy profile
$\beta(r)$ is utilized, which characterizes the velocity distribution’s
deviation from isotropy. This profile $\beta(r)$ is expressed as $1-(\langle
v_{\theta}^{2}\rangle+\langle v_{\phi}^{2}\rangle)/2\langle v_{r}^{2}\rangle$,
where $v_{\theta}$, $v_{\phi}$, and $v_{r}$ are the longitudinal, azimuthal,
and radial velocity components of galaxies, respectively (e.g. Diaferio
(1999)). The velocity anisotropy is generally an unknown variable in
kinematics analyses, and it should be correctly modeled along with the other
dynamical quantities. Parametric methods (e.g. Mamon  _et al._ (2013); Read
_et al._ (2021); Biviano  _et al._ (2023)) assume specific - physically
motivated - profiles which are fitted to the observed velocity and position
fields. Kinematics determinations based on the Jeans’ equation lead to a
degeneracy between $\beta(r)$ and the total cluster mass (the so-called mass-
anisotropy degeneracy, MAD); if additional information on the mass profiles is
provided (such as from X-ray or Weak Lensing analyses), non-parametric
reconstructions can be performed by inversion of the data (e.g. Binney and
Mamon (1982), see also Mamon  _et al._ (2019) and references therein).

The gravitational potential profile is related to the caustic amplitude
through the function $g(\beta)$, given by:

$g(\beta)=\frac{3-2\beta(r)}{1-\beta(r)}.$ (1)

As such, after estimating the caustic amplitude, $\beta(r)$ becomes the sole
unknown factor in estimating the gravitational potential. It’s worth
emphasizing that the caustic technique doesn’t depend on assumptions regarding
dynamical equilibrium, the shape of $g(\beta)$, or the gravitational potential
profile $\Phi(r)$ individually. Instead, it quantifies the combined effect of
$g(\beta)$ and $\Phi(r)$ in terms of the caustic amplitude
$\mathcal{A}^{2}(r)$,

$-2\Phi(r)=\mathcal{A}^{2}(r)g(\beta).$ (2)

This equation implies that the caustic technique can estimate a combination of
the gravitational potential profile and the velocity anisotropy parameter
$\beta$. This provides valuable insights into the dynamical properties of
spherical systems, particularly in the context of galaxy clusters and their
mass profiles.
The equation for the caustic mass profile of a spherical system can be
expressed as:

$GM(<r)=\int_{0}^{r}\mathcal{A}^{2}(r)\mathcal{F}_{\beta}(r)dr$ (3)

where $\mathcal{F}_{\beta}(r)=\mathcal{F}(r)g(\beta)$ and

$\mathcal{F}(r)=-2\pi G\frac{\rho(r)r^{2}}{\Phi(r)}.$ (4)

In this context, $\rho(r)$ is the density profile of the spherical system
(which we assume to be the NFW profile) and $\Phi(r)$ stands for the
gravitational potential profile.
Equation 3 relates the mass profile to the density profile of a spherical
system and a profile cannot be inferred without knowing the other. In
hierarchical clustering scenarios, the function $\mathcal{F}(r)$ doesn’t
exhibit strong variations with respect to radial distance $r$ Diaferio and
Geller (1997); Diaferio (1999). Similarly, $\mathcal{F}_{\beta}(r)$ also
changes slowly with $r$ if the velocity anisotropy parameter $\beta$ is
governed by a slowly varying function $g(\beta)$.
The practical utility of this technique becomes evident when applied to
individual clusters. It typically provides accurate escape velocity and mass
profiles, though there may be occasional deviations from the actual profile.
This method holds significant relevance as an alternative to gravitational
lensing for measuring mass in a cluster’s outer regions Diaferio  _et al._
(2005). Unlike lensing, it can be applied to clusters at any redshift,
provided there are a sufficient number of galaxies for a proper redshift
diagram analysis.
An important milestone in the application of the caustic technique was the
work of Geller, Diaferio, and Kurtz in 1999 Geller  _et al._ (1999). They
employed this method to measure the mass profile of the Coma cluster,
extending their analysis to an impressive $10h^{-1}$Mpc from the cluster
center. Their findings demonstrated that the Navarro, Frenk, & White (NFW)
profile provides an excellent fit to the cluster density profile at these
extensive radii, thereby challenging the viability of the isothermal sphere as
a cluster model. Subsequently, Biviano and Girardi (2003) Biviano and Girardi
(2003) applied the caustic technique to a composite cluster, stacking data
from 43 clusters in the Two Degree Galaxy Redshift Survey (2dGFRS, Colless et
al. 2001). Their results also aligned with earlier lensing and X-ray analyses,
affirming the method’s robustness and accuracy in estimating cluster mass
profiles.

##  III Modified gravity and Screening Mechanisms

In this section, we briefly describe the modified gravity scenarios assessed
in the current work and the underlying screening mechanisms.

Chameleon Screening: The chameleon model Khoury and Weltman (2004a);
Zaregonbadi  _et al._ (2022); Tsujikawa  _et al._ (2009); Kraiselburd  _et
al._ (2018); Khoury and Weltman (2004b) modifies gravity by introducing a
scalar field non-minimally coupled with the matter components and gives rise
to a fifth force that can be of the same order as the standard gravitational
force. The chameleon mechanism is a two-parameter model. The chameleon model
parameter $\beta$ determines the strength of the fifth force when it is not
screened. The second chameleon parameter, $\phi_{\infty}$, controls the
effectiveness of the screening mechanism, describing the transition from the
inner region of a cluster where gravity may be Newtonian to the outer region
where the fifth force contributes. The critical radius, where the transition
occurs, is determined by both $\phi_{\infty}$ and $\beta$. In the absence of
environmental effects, $\phi_{\infty}$ can be regarded to be the cosmological
background value of the chameleon field.

The chameleon mechanism operates whenever a scalar field couples to matter in
such a way that its effective mass depends on the local matter density. The
scalar-mediated force between matter particles can be of gravitational
strength, but its range is a decreasing function of ambient matter density,
thereby avoiding detection in regions of high density. Deep in space, where
the mass density is low, the scalar is light and mediates a fifth force of
gravitational strength, but near the Earth, where experiments are performed,
and where the local density is high, it acquires a large mass, making its
effects short range and hence, unobservable. This is achieved with a canonical
scalar field with suitable self-interaction potential $V(\phi)$, which is a
decreasing function of $\phi$ Khoury and Weltman (2004b). The theory (in the
weak-field limit and for non-relativistic matter) is given by Terukina  _et
al._ (2014); Khoury (2013):

$L_{\mathrm{chameleon}}=-\frac{1}{2}(\partial\phi)^{2}-V(\phi)-\frac{g\phi}{M_{\mathrm{Pl}}}\rho_{\mathrm{m}}$
(5)

The dimensionless coupling parameter $g$ is assumed to be $O(1)$,
corresponding to the gravitational strength coupling. In regions of high
density, the mass of the chameleon field increases. The range of interaction
decreases, thus screening the effect of the fifth force and recovering GR. In
low-density environments, the fifth force is unscreened and the effects of the
chameleon field can be observed. The modified gravitational potential under
this mechanism is given by Terukina  _et al._ (2014); Wilcox  _et al._ (2015):

$\frac{{\rm d}{\Phi(r)}}{{\rm
d}{r}}=\frac{G_{\mathrm{N}}{M}(r)}{r^{2}}+\beta\frac{{\rm d}{\phi}}{{\rm
d}{r}}$ (6)

where $\phi(r)$ is the chameleon field.
The chameleon field mediates a long-range fifth force when the matter density
is still large compared to the background, and the scalar field has not
settled in the minimum of the effective potential. The chameleon field becomes
effective beyond a critical radius, $r_{c}$, below which, the field is
completely screened. This radius is determined by Terukina  _et al._ (2014):

$1+\frac{r_{c}}{r_{s}}=\frac{\beta\rho_{s}r_{s}^{2}}{M_{\mathrm{Pl}}\phi_{\infty}}.$
(7)

The chameleon field Terukina  _et al._ (2014) in the two limits is given by:

$\phi(r)=\begin{cases}\phi_{s}[r/r_{s}(1+r/r_{s})^{2}]\equiv\phi_{int}(\approx
0)&(r<r_{c})\\\
-\frac{\beta\rho_{s}r_{s}^{2}}{M_{\mathrm{Pl}}}\frac{\mathrm{ln}(1+r/r_{s})}{r/r_{s}}-\frac{C}{r/r_{s}}+\phi_{\infty}\equiv\phi_{out}&(r>r_{c})\end{cases}$
(8)

where

$C\approx-\frac{\beta\rho_{s}r_{s}^{2}}{M_{\mathrm{Pl}}}\mathrm{ln}(1+r_{c}/r_{s})+\phi_{\infty}r_{c}/r_{s},$
(9)
$\phi_{\infty}-\frac{\beta\rho_{s}r_{s}^{2}}{M_{\mathrm{Pl}}}(1+r_{c}/r_{s})^{-1}\approx
0.$ (10)

Vainshtein Screening: In the case of Vainshtein screening the additional
degrees of freedom are screened through a non-linear mechanism. In Vainshtein
screening the gravitational potentials are modified inside the matter sources.
These modifications are screened outside the sources and GR is recovered in
low-density environments, as required by tests of gravity. GR is not recovered
everywhere within the Vainshtein radius as the screening breaks down there.
The Vainshtein radius $r_{V}$ is given by the curvature of the object Paillas
_et al._ (2019); Brando  _et al._ (2023).

The modified gravitational potential is now given as Cardone  _et al._ (2020);
Haridasu  _et al._ (2023); Crisostomi and Koyama (2018); Dima and Vernizzi
(2018); Langlois  _et al._ (2018); Bartolo  _et al._ (2018),

$\frac{{\rm d}\Phi(r)}{{\rm
d}r}=\frac{G_{\mathrm{N}}^{\mathrm{eff}}{M}(r)}{r^{2}}+\Xi_{1}G_{\mathrm{N}}^{\mathrm{eff}}{M}^{\prime\prime}(r)$
(11)

where $G_{\mathrm{N}}^{\mathrm{eff}}$ is the modified Newton’s constant
defined as $G_{\mathrm{N}}^{\mathrm{eff}}=\tilde{\gamma}_{\mathrm{N}}\times
G_{\rm N}$ and ′ represents the derivative w.r.t $r$. The effective parameter
of the theory $\Xi_{1}$ is related to the physical parameters as Langlois  _et
al._ (2017); Cardone  _et al._ (2020),

$\Xi_{1}=-\frac{(\alpha_{1}+\beta_{H})^{2}}{2(\alpha_{H}+2\beta_{1})}$ (12)

In Haridasu  _et al._ (2023), we set $\tilde{\gamma}_{\rm N}=1$ when
presenting the final results, as this parameter will be completely degenerate
with normalization of mass term ${M}(r)$. Hence, the assessment was limited
only to the parameter $\Xi_{1}$ as the physical parameters can only be
disentangled when the weak lensing potential is utilized complementary to the
hydrostatic potential Cardone  _et al._ (2020); Haridasu  _et al._ (2023), in
the case of Vainshtein screening, the gravitational potential and the weak
lensing potential are different. In the current work, although a joint
analysis with different mass-estimating techniques is utilized, the potential
is the same, unable to break additional degeneracies. However, we anticipate
that the improvement in the constraints of mass should help with better
constraining the $\Xi_{1}$ parameter. Although the $\tilde{\gamma}_{\rm N}$ is
unconstrained we leave it as a free parameter when performing the joint
importance sampling analysis, to assess the improvement in the constraints.
Note that both for the Chameleon and the Vainshtein screening mechanisms we
are solely interested in the modifications to the gravitational potential
given by Equations 11 and 6 to perform the comparisons with the caustic
surface estimated through the phase space of galaxies.

##  IV Data

In this work, we analyze X-ray and spectroscopic data of two massive galaxy
clusters, A2029 and A2142 at redshifts, $z=0.0784$ and $z=0.08982$,
respectively.

Hydrostatic data: These clusters have been amply studied, establishing their
hydrostatic masses and other thermodynamic properties, within the X-COP
compilation presented in Ettori  _et al._ (2019a); Eckert  _et al._ (2019);
Ghirardini  _et al._ (2019), which have been earlier utilized to test modified
gravity scenarios Ettori  _et al._ (2019a); Haridasu  _et al._ (2021, 2023);
Gandolfi  _et al._ (2023); Benetti  _et al._ (2023). Note that these two
clusters are relaxed Ettori  _et al._ (2019a) and have no non-thermal pressure
support at large radii Eckert  _et al._ (2019).

Caustic data: Meanwhile, the necessary data required for estimating caustics
are sourced from two references, namely Sohn  _et al._ (2017) for A2029 and
Liu  _et al._ (2018) for A2142. These datasets, detailed in Figure 2 and
Figure 3, consist of galaxy right ascension, declination, and redshift
information within the A2029 and A2142 clusters. This data forms the
foundation for constructing the velocity profiles of these clusters, as
visually represented in Figure 2 and Figure 3. For A2142 we have 2239 galaxies
and for A2029 we have 982 galaxies. The membership selection of these galaxies
for the formation of the caustics is done using the method detailed in Gifford
_et al._ (2013).

In effect, we choose to work with these two clusters as they have very good
hydrostatic and caustic data simultaneously. To calculate the gravitational
potentials associated with GR, CS, and VS models, as well as the resulting
profile $\mathcal{F}_{\beta}(r)$ from the hydrostatic data, we employ the
Markov Chain Monte Carlo (MCMC) chains, as utilized to constrain the models in
Boumechta  _et al._ (2023); Haridasu  _et al._ (2023, 2021). These chains
contain all the essential parameters such as $M_{200}^{\mathrm{HS}}$,
$R_{200}^{\mathrm{HS}}$, $c_{200}^{\mathrm{HS}}$, and the modified gravity
parameters.

##  V Method

Under the assumption that the cluster is mostly dominated by dark matter,
specifically in the range beyond $\sim 30\,\mathrm{kpc}$ which we consider, we
can model the mass density using the NFW profile Schaller  _et al._ (2015);
Hogan  _et al._ (2017); Sartoris  _et al._ (2020), which is given as Navarro
_et al._ (1996),

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

where $\rho_{s}$ is the characteristic density, $r_{s}$ is the characteristic
radius, where the logarithmic slope $s=d\mathrm{ln}\rho/d\mathrm{ln}r$ takes
the isothermal value s=-2. The corresponding mass profile Cardone  _et al._
(2020); Haridasu  _et al._ (2023) in terms of the $M_{\Delta}$ and
concentration $c_{\Delta}$ is given as,

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

where $x=r/R_{\Delta}$, $c_{\Delta}=R_{\Delta}/r_{s}$ and

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

We appropriately assume the value of the $\Delta$, to be to the usual value
$\Delta=200$ when making comparisons. we utilize the NFW mass profile to
assess the masses in the GR scenario and also in the Vainshtein and Chameleon
screening scenarios, as we have earlier done in Haridasu  _et al._ (2021,
2023) and Boumechta  _et al._ (2023), respectively. As the caustic technique
to obtain the mass profile is essentially independent of such an assumption on
the mass profile, it could be an interesting exercise to assess the same for
different assumptions of the mass profile, which we leave for a future
discussion.

###  V.1 Computing $\mathcal{F}_{\beta}(r)$

In the computation of $\mathcal{F}_{\beta}(r)$, we rely on the expressions
detailed in Equation 4 and Equation 1. Here, $\rho(r)$, as defined in Equation
13, serves as a crucial component and $\Phi(r)$ represents the gravitational
potential. In the context of General Relativity (GR), $\Phi(r)$ corresponds to
the familiar Newtonian potential. In the Chameleon Scalar-Tensor (CS) and
Vector Scalar-Tensor (VS) models, $\Phi(r)$ assumes the forms outlined in
Equation 6 and Equation 11, respectively.
We perform the calculation for $\mathcal{F}_{\beta}(r)$ across 1000 profiles,
each corresponding to distinct parameter values drawn from the MCMC111We
utilize the emcee code Foreman-Mackey  _et al._ (2013) publicly available at
https://github.com/dfm/emcee. samples obtained after fitting the hydrostatic
data, see Boumechta  _et al._ (2023) for Chameleon, and Haridasu  _et al._
(2023) for Vainshtein screening. In the resulting plot depicted in Figure 1,
the solid lines represent the mean values extracted from these profiles, while
the shaded region encompasses the $95\%$ C.L. allowed regions of the profiles.

###  V.2 Measuring $\mathcal{A}(r)$

To determine the caustic amplitude, denoted as $\mathcal{A}^{2}(r)$, we employ
a methodology based on the velocity dispersion profile derived from
observational data. This technique, as described in Diaferio (1999), relies on
the velocity distribution of galaxy clusters. We begin by utilizing velocity
data obtained from the members of galaxy clusters, as provided in the datasets
Sohn  _et al._ (2017) and Liu  _et al._ (2018) for A2029 and A2142
respectively. Following the selection of velocities Gifford  _et al._
(2013)222We use part of the code provided in Gifford  _et al._ (2013) to
produce the caustic surfaces and extract the caustic amplitude for the mass
calculation., taking into account the full phase-space Wing and Blanton (2011)
for inclusion in the caustic diagram, we proceed to define the caustic
surfaces using the technique detailed in Diaferio (1999). Consider a
collection of $N$ galaxies, each characterized by coordinates
$\mathbf{x}=(r,v)$, where we conveniently rescale both $r$ and $v$. We adopt
an adaptive kernel method to estimate the density distribution of these
galaxies within the redshift diagram.

This density is then expressed as

$f_{q}(\mathbf{x})=\frac{1}{N}\sum_{i=0}^{N}\frac{1}{h_{i}^{2}}K\left(\frac{\mathbf{x}-\mathbf{x}_{i}}{h_{i}}\right).$
(16)

And the kernel function $K(t)$ is defined as

$K(t)=\begin{cases}4\pi^{-1}(1-t^{2})^{3}&t<1\\\
1&\mathrm{otherwise}.\end{cases}$ (17)

The local smoothing length, denoted as $h_{i}$, depends on both the local
density and an optical smoothing length $h_{\mathrm{opt}}$, which is given by:

$h_{\mathrm{opt}}=\frac{3.12}{N^{1/6}}\left(\frac{\sigma_{r}^{2}+\sigma_{v}^{2}}{2}\right)^{1/2}.$
(18)

Here, $\sigma_{r}$ and $\sigma_{v}$ represent the marginal standard deviations
of the galaxy coordinates. Additionally, a local smoothing factor
$\lambda_{i}$ is introduced as
$\lambda_{i}=[\gamma/f_{1}(\mathbf{x}{i})]^{1/2}$, where $f_{1}$ is calculated
from the density distribution. Notably, we set $h_{c}=\lambda_{i}=1$ for any
$i$, and $\log\gamma$ is computed as $\sum_{i}\log[f_{1}(\mathbf{x}_{i})]/N$.

The amplitude $\mathcal{A}^{2}(r)$ at a fixed radial distance $r$ is then
determined by solving the equation $f_{q}(r,v)=\kappa$. More specifically, the
first upper and lower solutions, denoted as $v_{\mathrm{u}}$ and
$v_{\mathrm{d}}$, are identified away from the maximum of $f_{q}(r,v)$,
closest to $v=0$. The amplitude $\mathcal{A}(r)$ is computed as
$\min\\{|v_{\mathrm{u}}|,|v_{\mathrm{d}}|\\}$. It’s worth noting that the
prescription $\mathcal{A}(r)=\min{|v_{u}|,|v_{d}|}$ is equivalent to
$\mathcal{A}(r)=(v_{\mathrm{u}}-v_{\mathrm{d}})/2$ in the case of an isolated
spherically symmetric system.

While we obtain $f_{q}(r,v)$ uniquely except for the choice of $q$ 333The
value of $q$ could be assumed between 5-10, and we validate that it does not
affect the estimation of the caustic profile, within the estimated error. In
the current work, we utilize $q=10$, finding no major difference with the
assumed value. , it’s important to highlight that there exists an infinite
range of thresholds $\kappa$ that can be employed to determine
$\mathcal{A}(r)$. In practice, it’s reasonable to assume that, particularly in
the central region, the cluster has reached a state of virial stability. Thus,
within this central region, the equation $\langle
v^{2}_{\mathrm{esc}}\rangle_{R}=\langle v^{2}\rangle_{R}$ holds. Here, the
angular brackets indicate an average computed over the entire sphere of radius
$R$, considering velocities in three dimensions. In our specific dataset, for
example, from Sohn  _et al._ (2017), we only have access to one-dimensional
velocity information. Therefore, the methodology also assumes that if the
velocity field exhibits approximate isotropy in the central region, our
expression remains valid. Consequently, we have $\langle
v^{2}\rangle_{\kappa,R}=\int_{0}^{R}\mathcal{A}^{2}(r)\Phi(r)dr/\int_{0}^{R}\Phi(r)dr$,
where $\Phi(r)$ represents the integral of $f_{q}(r,v)$. It’s important to
emphasize that $\langle v^{2}_{\mathrm{esc}}\rangle_{\kappa,R}$ is the only
quantity that depends on the chosen threshold $\kappa$.

To select an appropriate $\kappa$, we minimize the function defined as:

$S(\kappa,R)=|\langle v_{\mathrm{esc}}^{2}\rangle_{\kappa,R}-4\langle
v^{2}\rangle_{R}|^{2}.$ (19)

Here, we define the radial distance $R$ as the mean projected distance of
cluster members from the cluster center, and $\langle v^{2}\rangle_{R}$
represents the one-dimensional velocity dispersion of the cluster members.
This comprehensive methodology enables us to robustly estimate the caustic
amplitude and, by extension, gain insights into the cluster’s dynamical
properties and mass distribution.

The uncertainty in the measured value of $\mathcal{A}(r)$ depends on the
number of galaxies contributing to the determination of $\mathcal{A}(r)$ Serra
_et al._ (2010). Therefore, we define the relative error,

$\delta\mathcal{A}(r)/\mathcal{A}(r)=\kappa/\mathrm{max}\\{f_{q}(r,v)\\}$ (20)

where the maximum value of $f_{q}(r,v)$ is along the $v$-axis at fixed $r$.
The resulting error on the cumulative mass profile is

$\delta
M_{i}=\sum_{j=1}^{i}|2m_{j}\delta\mathcal{A}(r_{j})/\mathcal{A}(r_{j})|$ (21)

where $m_{j}$ is the mass of the shell $\\{r_{j-1},r_{j}\\}$ and $i$ is the
index of the most external shell.

Figure 1:  We show the profiles of $\mathcal{F}_{\beta}(r)$ constructed
utilizing the constraints on the gravitational potential obtained by fitting
the hydrostatic data. In the top we show the profiles for the cluster A2029,
for GR (green), Chameleon screening (blue), and Vainshtein screening (red).
Similarly, in bottom panel we show the same for the cluster A2142.

###  V.3 Comparing $\mathcal{A}^{2}(r)$

As a preliminary assessment, we utilize Equation 2, to perform a comparison
between the gravitational potential on the L.H.S estimated using the
hydrostatic equilibrium and the caustic surface on R.H.S, using the phase
space data of the galaxies. The fundamental definition of $\mathcal{A}^{2}(r)$
is articulated in Equation 2. In our investigation, we undertake a comparative
study involving caustic profiles derived from empirical observations and
theoretical projections of $\mathcal{A}^{2}(r)$. The latter is determined
through both Equation 2 and the inversion of Equation 3. While Equation 2
yields a prediction for $\mathcal{A}^{2}(r)$ independent of
$\mathcal{F}_{\beta}(r)$, Equation 3 provides predictions contingent upon
$\mathcal{F}_{\beta}(r)$. We meticulously scrutinize how $\mathcal{A}^{2}(r)$
behaves under constant and varying $\mathcal{F}_{\beta}(r)$ scenarios.
Additionally, these equations allow us to quantitatively assess the
fluctuations in $\mathcal{A}^{2}(r)$ within the context of modified
gravitational potentials while assuming an NFW mass profile, as defined in
Equation 14. This approach empowers us to estimate the requisite mass
adjustments to align a modified gravitational potential with the observed
$\mathcal{A}^{2}(r)$. In essence, this methodology provides an indirect route
for evaluating the impact of modified gravity on mass bias, an essential
aspect of our study, serving as both a justification and validation of our
findings. A direct method for obtaining caustic mass involves the
straightforward application of Equation 3. This approach constitutes the core
of our caustic mass determination.

The computation of caustic profiles and the assessment of caustic amplitude
primarily rely on the velocities and velocity dispersion of galaxy clusters,
as outlined by Gifford  _et al._ (2013). Detailed insights into the
methodology for obtaining caustic amplitude, denoted as $\mathcal{A}^{2}(r)$,
can be found in Diaferio (1999) and Gifford  _et al._ (2013). In the existing
literature, $\mathcal{F}_{\beta}(r)$ is typically approximated as either 0.5
Diaferio (1999) or 0.65 Gifford  _et al._ (2013). In contrast, we implement a
more comprehensive approach by utilizing the complete expression for
$\mathcal{F}_{\beta}(r)$ involving Equation 4 and Equation 1 apriori derived
independently from the hydrostatic data, which in turn is the important aspect
of our joint analysis.

In GR scenarios, the gravitational potential generated by the cluster,
$\Phi(r)$, is commonly known to be Diaferio (1999),

$\Phi(r)=-\frac{G{M}(<r)}{r}-4\pi G\int^{\infty}_{r}\rho(x)x{\rm d}x$ (22)

and $\rho(r)$ is the NFW density profile given in Equation 13. It is
straightforward to evaluate $\mathcal{F}(r)$ in this scenario:

$\mathcal{F}_{\mathrm{NFW}}(r)=\frac{r^{2}}{2(r+r_{s})^{2}}\frac{1}{\mathrm{ln}(1+r/r_{s})}$
(23)

where $r_{s}$ is the scale radius Diaferio (2009). Furthermore, we extend our
analysis to modified gravity scenarios by evaluating $\mathcal{F}_{\beta}(r)$
using our modified potentials, Equation 6 and Equation 11. This multifaceted
approach equips us to comprehensively explore the intricate interplay between
mass predictions, gravitational potentials, and modified gravity scenarios. In
Figure 2 and Figure 3 we present a comparative analysis of the caustic
profiles derived from the two methodologies. In each plot, the blue line
labels the caustic profile obtained with the caustic technique, while the
green line, labeled as Hydrostatic-GR/CS/VS, shows the caustic profile
independently predicted by the hydrostatic equilibrium analysis for the three
models, respectively.

Remarkably, the caustic profiles obtained through these two methods exhibit a
commendable level of agreement. The caustic profiles estimated using a
constant $\mathcal{F}_{\beta}(r)=0.65$ also agree reasonably well except for
the innermost regions. Similarly, an $\mathcal{F}_{\beta}(r)=0.5$
underestimates the caustic profile, and consequently, we can anticipate that
the mass estimated assuming $\mathcal{F}_{\beta}(r)=0.5$ would also be an
underestimation. In the current analysis, alongside our method to utilize the
hydrostatic data we also assume $\mathcal{F}_{\beta}(r)=0.5$ or
$\mathcal{F}_{\beta}(r)=0.65$ to illustrate the differences and further
discuss the implications, especially in the modified gravity scenarios. While
an analysis to anticipate a $\mathcal{F}_{\beta}(r)$ in the modified gravity
scenarios is not present, a recent analysis Pizzardo  _et al._ (2023a) in the
case of GR find $\mathcal{F}_{\beta}(r)=0.4\pm 0.08$ using simulations
Springel  _et al._ (2018), which is in agreement with the $0.5$ value. To
further illustrate potential deviations in the caustic profiles and - by
extension - the mass profiles induced in the non-GR case, in the central and
left plot of each figure of Figures 2 and 3 we additionally show caustic
profiles for a few specific values of the modified gravity parameters in both
Vainshtein and Chameleon screening.

Figure 2: Caustic Profiles for the cluster A2029. We show the caustic profile
and the corresponding uncertainty in shaded blue regions, alongside the phase
space distribution of the galaxies. Left: We compare the caustic profiles
estimated using the hydrostatic data in the case of GR, shown in green. The
caustic profile estimated using the hydrostatic data and a constant
$\mathcal{F}_{\beta}(r)=0.65$ is shown in red Center: Same as the Left panel
comparing against the caustic surface estimated using the hydrostatic data in
the case of chameleon screening. Right: Same as the other two panels,
comparing the caustic surface in the case of Vainshtein screening. For the
Chameleon and Vainstein screening cases, we show a few expected caustic
surfaces for various values of assumed values of the modified gravity
parameters.

Figure 3: Same as Figure 2, but for the cluster A2142. Note the change in the
radial range with respect to cluster A2029.

##  VI Results and Discussion

We begin by presenting our results for the construction of the caustic
profiles for the two clusters considered in our analysis. In Figures 2 and 3,
we show the caustic profiles of the clusters A2029 and A2142, respectively.
Alongside the reconstructed caustic surface (blue) we also show the caustic
surface as anticipated by the hydrostatic technique assuming an NFW mass
profile (green), which essentially is a first step validation to further
proceed with the joint analysis. In general, we notice a very good agreement
between the two estimates. In the Center and the Right panels we show the same
comparison against the Chameleon and Vainshtein screening, respectively. In
the case of Chameleon screening, especially for cluster A2029 one can notice a
sharp deviation of the caustic surface predicted by hydrostatic equilibrium at
the screening radius. It is worth pointing out that this radial range is
fairly outside the range covered by hydrostatic data and cannot be
constrained, however when compared with the caustic surface one can easily
rule out the certain parameter space that provides such a large deviation.
Essentially, this demonstrates the added advantage of the joint analysis
including the caustic technique that can probe the outskirts of the galaxy
cluster. In the case of the Vainshtein screening, however, the modifications
are not limited to the radial range larger than the screening radius, as in
this case, the fifth force penetrates the inner regions of the galaxy
clusters. This can be easily understood by contrasting the modifications of
the gravitational potential in these two modified gravity scenarios (see
Equations 11 and 6). Needless to say, we begin by contrasting our caustic mass
estimates, to those present in literature to validate our procedure assuming a
constant $\mathcal{F}_{\beta}(r)$.

###  VI.1 Joint constraints on Mass

Firstly, we assess the constraint on mass in the standard GR case taking into
account the appropriate $\mathcal{F}_{\beta}(r)$, well-constrained using the
hydrostatic data. In Figure 4, we show the marginalized posteriors for the
hydrostatic mass compared against the caustic mass when assuming a constant
$\mathcal{F}_{\beta}(r)$, which are also reported in Table 1. As it can be
seen from the contours, the caustic mass evaluated at
$R_{200}^{\rm{Cau}}$444The $R_{200}^{\rm{Cau}}$ is simply evaluated as the
radius within which the average density is 200 times the critical density,
after the mass profile is reconstructed. is mildly anti-correlated with the
$M_{200}^{\rm{HS}}$, while the mass estimated at $R_{200}^{\rm{HS}}$ from the
caustic mass profile shows a positive correlation, which is an expected
behavior. Once we perform the importance sampling all the mass estimates are
positively correlated and heavily affected by the $\mathcal{F}_{\beta}(r)$
derived using the hydrostatic data.

We find that the mass bias earlier quoted in the comparison of the hydrostatic
mass and the caustic mass is completely alleviated in our joint analysis. For
instance, Ettori  _et al._ (2019b) quote a $6.5\sigma$ and $3.9\sigma$
discrepancy between the two masses, for the clusters A2029 and A2142,
respectively. These clusters are among the high significance of deviation
analyses within the X-COP compilation. However, this should not come as a
surprise as it is well known that the global assumption of
$\mathcal{F}_{\beta}(r)=0.5$ Diaferio  _et al._ (2005) or
$\mathcal{F}_{\beta}(r)=0.65$ Gifford  _et al._ (2013), can lead to such a
bias. As can be seen in Table 1, yielding lower mass values, the assumption of
$\mathcal{F}_{\beta}(r)=0.5$ leads to a larger bias w.r.t. the
$\mathcal{F}_{\beta}(r)=0.65$ case. When taking the appropriate hydrostatic
$\mathcal{F}_{\beta}(r)$ we find the mass bias to be no more than $\sim
0.03-0.1$, at a significance of $\lesssim 1\sigma$ for both the clusters (see
the last column of Table 1). For the sake of comparison, we also show the mass
reconstructed using constant $\mathcal{F}_{\beta}(r)$ and our method in fig.
9. Incidentally, the assumption of $\mathcal{F}_{\beta}(r)=0.65$, yields a
very good agreement between the independent hydrostatic and caustic mass
estimates. In the case of A2142 the assumption of
$\mathcal{F}_{\beta}(r)=0.65$ is in fact in slightly better agreement than our
analysis using the hydrostatic data.

On the other hand, we also find that the caustic mass is mildly larger than
the hydrostatic masses. For ease of comparison, we present the posteriors of
the masses within the GR case in Figure 5. Given the very good agreement
between the hydrostatic mass and the caustic mass obtained within our
formalism, we now perform an importance sampling analysis to obtain the joint
constraint on the mass of the cluster. The corresponding values of the
constraints on mass are reported in the penultimate column of Table 1 and
shown as red posteriors in Figure 5. Our joint analysis formalism allows us to
obtain tighter constraints on the mass of the clusters. To the best of our
knowledge, this is the first combined analysis between these two mass
estimation techniques, forward modeling of the hydrostatic data and caustic
technique, to obtain a joint constraint that now allows us to extend the
formalism to alternate theories of gravity, which indeed is the primary goal
of the current investigation.

Table 1: Constraints on the mass in the case of GR, shown as $68\%$ C.L. limits. In the second and third columns, we show the constraints obtained using the hydrostatic data only. In columns five and six we quote the caustic mass estimates and in column seven we show the joint mass estimate. In the last column, we show the mass bias between the hydrostatic and caustic techniques. Column 4 shows the corresponding $\mathcal{F}_{\beta}(r)$ assumption used in the analysis to estimate the caustic mass.  Cluster | $R_{200}^{\rm{HS}}$ | $M_{200}^{\rm{HS}}$ | $\mathcal{F}_{\beta}(r)$ | $R_{200}^{\rm{Cau}}$ | $M_{200}^{\rm{Cau}}$ | $M_{200}^{\rm{Joint}}$ |  $M_{200}^{\rm{HS}}$/$M_{200}^{\rm{Cau}}$
---|---|---|---|---|---|---|---
| [Mpc] | $[10^{14}\,M_{\odot}]$ |  | [Mpc] | $[10^{14}\,M_{\odot}]$ | $[10^{14}\,M_{\odot}]$ |
A2029 | $2.08_{-0.03}^{+0.03}$ | $11.04_{-0.47}^{+0.50}$ | HS | $2.11_{-0.03}^{+0.03}$ | $11.43_{-0.54}^{+0.55}$ | $11.21_{-0.34}^{+0.33}$ | $0.97_{-0.07}^{+0.07}$
0.65 | $2.02_{-0.03}^{+0.03}$ | $10.25_{-0.49}^{+0.51}$ | $10.79_{-0.38}^{+0.39}$ | $1.08_{-0.07}^{+0.07}$
0.50 | $1.86_{-0.03}^{+0.03}$ | $7.89_{-0.38}^{+0.39}$ | $9.76_{-0.11}^{+0.05}$ | $1.40_{-0.09}^{+0.09}$
A2142 | $2.15_{-0.03}^{+0.03}$ | $12.39_{-0.52}^{+0.57}$ | HS | $2.23_{-0.03}^{+0.06}$ | $13.70_{-0.95}^{+1.06}$ | $12.75_{-0.43}^{+0.43}$ | $0.90_{-0.08}^{+0.09}$
0.65 | $2.20_{-0.04}^{+0.04}$ | $13.20_{-0.79}^{+0.84}$ | $12.63_{-0.44}^{+0.47}$ | $0.94_{-0.07}^{+0.07}$
0.50 | $2.02_{-0.04}^{+0.04}$ | $10.16_{-0.61}^{+0.65}$ | $11.62_{-0.40}^{+0.41}$ | $1.22_{-0.09}^{+0.10}$

Figure 4: Left: We show the $68\%$ and $95\%$, confidence levels of the
$M_{200}$ values, estimated using the hydrostatic data, assuming the constant
$\mathcal{F}_{\beta}(r)=0.65$ cases. The unfilled contours represent the final
importance-sampled posteriors corresponding to the penultimate column of Table
1. Figure 5: Comparison of the posteriors on the mass estimates from different
techniques and the joint analysis, for the clusters A2029 (top) and A2142
(bottom) panels.

###  VI.2 Constraints on Modified gravity

As mentioned in the previous section, finding very good agreement between the
hydrostatic and caustic masses inferred with our method in the GR case, now we
extend to the modified gravity case simultaneously assessing both mass bias
and the constraints on the additional fifth force parameters. In Tables 2 and
3, we compare the constraints on mass and modified gravity parameters. We show
a simple comparison of the constraints when assuming
$\mathcal{F}_{\beta}(r)=0.5$, to highlight the scenario wherein the caustic
mass is inferred independent of the hydrostatic mass, and then a joint
analysis is performed. This is a valid scenario given that the masses in the
modified gravity models are less constrained and allow posteriors that can be
consistent with $\mathcal{F}_{\beta}(r)=0.5$ and hence can provide varied
constraints for the modified gravity parameters (see Appendix C). Note that
this comparison should be taken at face value only for illustrative purposes.
In Figures 6 and 7 we show the confidence level of the posteriors for the
same, for the modified gravity parameters.

Table 2: CS Masses Cluster | $R_{200}^{\rm{HS}}$ | $M_{200}^{\rm{HS}}$ | $\mathcal{F}_{\beta}(r)$ | $R_{200}^{\rm{Cau}}$ | $M_{200}^{\rm{Cau}}$ | $M_{200}^{\rm{Joint}}$ | $\phi_{\infty,2}$ | $\beta_{2}$ |  $M_{200}^{\rm{HS}}$/$M_{200}^{\rm{Cau}}$
---|---|---|---|---|---|---|---|---|---
| [Mpc] | $[10^{14}\,M_{\odot}]$ |  | [Mpc] | $[10^{14}\,M_{\odot}]$ | $[10^{14}\,M_{\odot}]$ | $\beta=\sqrt{\frac{1}{6}}$ |  |
A2029 | $2.04_{-0.13}^{+0.04}$ | $10.41_{-1.91}^{+0.62}$ | HS | $2.08_{-0.10}^{+0.05}$ | $11.02_{-1.57}^{+0.75}$ | $10.84_{-0.58}^{+0.39}$ | $0.24_{-0.18}^{+0.29}$ | – | $0.94_{-0.12}^{+0.13}$
0.50 | $1.86_{-0.03}^{+0.04}$ | $7.83_{-0.37}^{+0.38}$ | $7.87_{-0.46}^{+0.42}$ | $0.36_{-0.03}^{+0.12}$ | $0.28_{-0.03}^{+0.03}$ | $1.32_{-0.23}^{+0.12}$
A2142 | $2.14_{-0.13}^{+0.03}$ | $12.20_{-2.13}^{+0.53}$ | HS | $2.20_{-0.13}^{+0.07}$ | $13.19_{-2.21}^{+1.20}$ | $12.48_{-0.74}^{+0.44}$ | $0.28_{-0.21}^{+0.35}$ | – | $0.91_{-0.10}^{+0.15}$
0.50 | $2.02_{-0.04}^{+0.04}$ | $10.18_{-0.65}^{+0.59}$ | $10.28_{-0.64}^{+0.78}$ | $0.45_{-0.12}^{+0.36}$ | $0.22_{-0.04}^{+0.04}$ | $1.18_{-0.19}^{+0.11}$

Figure 6: We show the $68\%$ and $95\%$ C.L. contours in the case of the
Chameleon screening, for the clusters A2029 (left) and A2142 (right). The
filled contours show the constraints obtained using only the hydrostatic data
and the open contours show the joint analysis including the caustic technique.

Chameleon screening: As can be seen in the Figure 6, the joint analysis of the
hydrostatic and caustic techniques helps reduce the degenerate region between
mass ($M_{500}$) and the coupling constant ($\beta_{2}$) parameter, which we
earlier elaborated upon in Boumechta  _et al._ (2023). Traditionally, the
approach to alleviate this degeneracy and obtain constraints was to assume a
mass prior obtained using weak lensing techniques and in Boumechta  _et al._
(2023), we introduced an internal mass prior that a prior of the mass obtained
when restricting the analysis to $\beta_{2}\geq 0.5$. In the current work, we
utilize the Caustic technique to fulfill the role of a prior, however, the
formalism is not to utilize a mass prior but to perform an importance sampling
while simultaneously constraining the mass from the caustic technique. The
ability of the Caustic technique to improve the constraints heavily relies on
the mass estimates alone. In Table 2, we show the joint mass estimates, which
are once again tighter than either of the individual masses. In the last
column of Table 2, we show the mass bias between the hydrostatic and the
caustic techniques, where we find a very good consistency within $\lesssim
1\sigma$.

Needless to say, as can be seen in Figure 6, the constrain the
$\phi_{\infty,2}$ for the case of $\beta_{2}=\sqrt{1/6}$, which corresponds to
the specific case of $f(R)$ gravity are also largely reduced for both the
clusters. We infer a joint constraint of $\phi_{\infty,2}<0.18$ at $99\%$
C.L., using only the two clusters here. This is a very good constraint and is
comparable to the joint analysis present in Boumechta  _et al._ (2023), where
five clusters were utilized with Weak lensing priors, $\phi_{\infty,2}<0.13$
at the same confidence level. This improvement here is essentially because the
caustic mass we constraint taking into account the $\mathcal{F}_{\beta}(r)$
constrained using the hydrostatic data is much tighter than the WL masses
present in Herbonnet  _et al._ (2020) (see Table I in Boumechta  _et al._
(2023)), which were utilized in Boumechta  _et al._ (2023) (see Table II
therein).

Vainshtein screening: In the case of Vainshtein screening, we have earlier
established the constrained on the $\Xi$ parameter in Haridasu  _et al._
(2023). As can be seen in the third column of Table 3, the current clusters in
A2029, A2142 are very good examples for GR and mild deviation from GR case,
respectively. In this context, the joint analysis is very helpful in
constraining the modified gravity parameter $\Xi_{1}$. Firstly, as described
in Section III (see also Haridasu  _et al._ (2023)), within the Vainshtein
screening the gravitational potential and the weak lensing potential are
distinct a joint analysis would be needed to correctly constrain the running
of the gravitational constant $G_{\rm N}^{\rm eff}$. As the
$\mathcal{F}_{\beta}(r)$ necessary to constrain the caustic mass is the
gravitational potential as well, we do not immediately break the degeneracy
between the $\tilde{\gamma}_{\rm N}$ and $M_{200}^{\rm{HS}}$. We present the
constraints on the combination of these parameters in Table 3.

In the case of cluster A2029, we find a twice as tighter constraint for the
parameter $\Xi_{1}=0.0\pm 0.10$ with respect to the constraint of
$\Xi_{1}=-0.04^{+0.19}_{-0.12}$ obtained using only the hydrostatic data. This
constraint is obtained while having no bias between the estimates, with
$M_{200}^{\rm{HS}}/M_{200}^{\rm{Cau}}=0.98^{+0.15}_{-0.11}$, which shows a
complete consistency. On the other hand, cluster A2142 which shows a mild
$\sim 2\sigma$ deviation with $\Xi_{1}\sim-0.20^{+0.10}_{-0.08}$ using only
the hydrostatic data, is now more consistent with GR within $\sim 1\sigma$,
having $\Xi_{1}=-0.10^{+0.11}_{-0.06}$. Note that this is a shift in the
parameters with no major improvement in the relative constraint. This
improvement in the agreement with the GR is essentially because the masses
estimated using the hydrostatic method and the caustic method have a bias of
$M_{200}^{\rm{HS}}/M_{200}^{\rm{Cau}}=0.81^{+0.09}_{-0.07}$, which is about
the same significance $\geq 2\sigma$, as the earlier for modified gravity
using the hydrostatic data alone. Nevertheless, we perform the joint analysis
as the inconsistency is only of the order of $\sim 2\sigma$.

In this context, it is also interesting to note that for the cluster A2142, in
the case of $\mathcal{F}_{\beta}(r)=0.5$, we find that the mass bias is no
longer at a $\sim 2\sigma$ significance as earlier and is completely
consistent with unity being
$M_{200}^{\rm{HS}}/M_{200}^{\rm{Cau}}=1.05^{+0.13}_{-0.11}$. This is, however,
accompanied by the joint constraint of $\Xi_{1}=-0.21\pm 0.07$, which is a
$3\sigma$ detection of modified gravity. We find this to be an interesting
observation: if one demands the masses estimated using the two techniques to
be consistent as a prior before assessing the joint constraint on the modified
gravity parameters, this hypothetical case where the assumption of
$\mathcal{F}_{\beta}(r)=0.5$ is not justified seems to provide better
agreement and detection of deviation from GR. At face value, not having a mass
bias, in this case, could imply that the extreme bias quoted in Ettori  _et
al._ (2019b) could be completely alleviated by simply changing the gravity
model with the assumption of $\mathcal{F}_{\beta}(r)$. However, it is
important to assess the $\mathcal{F}_{\beta}(r)$ correctly, before addressing
the mass bias, as we have done consistently in our formalism even though the
data shows a mild bias. As can be seen in Table 3, this is not the case for
cluster A2029, which is completely consistent with GR.

Table 3: VS Masses Cluster | $R_{200}^{\mathrm{HS}}$ | $\Xi_{1}$ | $\tilde{\gamma}_{\rm N}\,M_{200}^{\mathrm{HS}}$ | $\mathcal{F}_{\beta}(r)$ | $R_{200}^{\mathrm{Cau}}$ | $\tilde{\gamma}_{\rm N}\,M_{200}^{\mathrm{Cau}}$ | $\tilde{\gamma}_{\rm N}\,M_{200}^{\mathrm{Joint}}$ | $\Xi_{1}$ |  $M_{200}^{\rm{HS}}$/$M_{200}^{\rm{Cau}}$
---|---|---|---|---|---|---|---|---|---
| [Mpc] |  |  |  | [Mpc] | $[10^{14}\,M_{\odot}]$ | $[10^{14}\,M_{\odot}]$ |  |
A2029 | $2.09_{-0.09}^{+0.11}$ | $-0.04_{-0.12}^{+0.19}$ | $11.16_{-1.40}^{+1.92}$ | $\mathrm{HS}$ | $2.10_{-0.04}^{+0.04}$ | $11.32_{-0.64}^{+0.68}$ | $11.15_{-1.01}^{+1.15}$ | $-0.00_{-0.09}^{+0.10}$ | $0.98_{-0.11}^{+0.15}$
0.50 | $1.86_{-0.03}^{+0.03}$ | $7.84_{-0.38}^{+0.39}$ | $8.78_{-0.44}^{+0.53}$ | $-0.18_{-0.06}^{+0.06}$ | $1.43_{-0.19}^{+0.25}$
A2142 | $2.04_{-0.06}^{+0.07}$ | $-0.20_{-0.08}^{+0.10}$ | $10.50_{-0.87}^{+1.04}$ | $\mathrm{HS}$ | $2.19_{-0.05}^{+0.05}$ | $13.01_{-0.90}^{+0.96}$ | $11.80_{-0.94}^{+1.05}$ | $-0.10_{-0.06}^{+0.11}$ | $0.81_{-0.07}^{+0.09}$
$0.50$ | $2.01_{-0.04}^{+0.04}$ | $10.04_{-0.62}^{+0.64}$ | $10.21_{-0.58}^{+0.60}$ | $-0.21_{-0.06}^{+0.07}$ | $1.05_{-0.11}^{+0.13}$

Figure 7: We show the $68\%$ and $95\%$ C.L. contours in the case of the
Vainshtein screening, for the clusters A2029 (left) and A2142 (right). The
filled contours show the constraints obtained using only the hydrostatic data
and the open contours show the joint analysis including the caustic technique.

##  VII Conclusion

Joint analysis of different mass estimating techniques of galaxy clusters has
been time and again utilized in constraining accurately the mass and
consequently testing modified gravity theories Terukina  _et al._ (2014);
Wilcox  _et al._ (2015); Sakstein (2015); Haridasu  _et al._ (2023); Pizzuti
_et al._ (2020, 2023). In the same spirit, here we have assessed the ability
of the caustic technique Diaferio and Geller (1997); Diaferio (2009); Serra
_et al._ (2010) to aid in testing modified gravity scenarios for the first
time, to the best of our knowledge.

The main implications and findings are as follows:

  * •

Firstly, we validate our formalism in GR, obtaining joint constraints on the
mass of the two clusters. Essentially this is aided by the fact that we see no
mass bias, which is quoted in Ettori  _et al._ (2019a). While we have already
confirmed and found that the mass bias is reduced when the hydrostatic
$\mathcal{F}_{\beta}(r)$ is taken into account, another possible way to
alleviate the mass bias would be by assuming a modified gravity scenario.

  * •

We have shown that the Caustic technique can robustly aid in testing the
gravity within galaxy clusters. We have tested for the cases of Vainshtein and
Chameleon screening. We also perform a hypothetical analysis, in case the
caustic mass is estimated independently, how it could lead to spurious
detection of modified gravity.

  * •

In the case of the Vainshtein screening we were able to obtain tighter
constraints on the free parameter with respect to the solely hydrostatic
analysis. In particular, for the cluster A2142, the X-ray data provided a
constraint below the expectation of GR, $\Xi_{1}\sim-0.2$, with $\sim 2\sigma$
significance, while the joint analysis reduces the tension to
$\Xi_{1}=-0.10^{+0.11}_{-0.06}$ at $\lesssim 1\sigma$. For cluster A2029, the
consistency with GR remains intact with an improvement on the bounds at the
same confidence level of a factor of $\sim 2$ in the joint analysis.

  * •

For the case of Chameleon screening we find that the caustic technique aids in
reducing the degeneracy between the mass estimates ($M_{500}$) and the
coupling constant ($\beta_{2}$), as can be seen in Figure 6. Given this we
estimate a competitive joint constraint of $\phi_{\infty,2}\lesssim 0.18$ at
$3\sigma$ C.L., using only the two clusters utilized in this work. This
further reinforces the utility of the caustic technique in constraining the
modified gravity parameters.

Although in the current work, we have limited ourselves to two clusters as
proof of concept to explore the possibility, the formalism can be easily
extended to several individual clusters stacked Gifford  _et al._ (2017);
Svensmark  _et al._ (2015); Pizzardo  _et al._ (2023a, b) for both the
hydrostatic data (see for example Wilcox  _et al._ (2015, 2016)) and caustic
phase space, which we intend to extend as one of the future directions. We
anticipate the current formalism to help perform a thorough investigation of
the modified gravity models utilizing the caustics as an independent
observable alongside the already well-explored weak lensing and the
hydrostatic datasets. we also intend to extend our formalism to other modified
gravity, non-standard dark matter scenarios Gandolfi  _et al._ (2023); Benetti
_et al._ (2023), with non-local effects on the dark matter mass profiles, in
future communications. Aside from these we also intend to perfrom a joint
analysis of the kinematics, dynamics, and weak lensing datasets exploring the
potential to constrain the modified gravity parameters and break degeneracies
within them.

## Acknowledgements

The authors are grateful to Antonaldo Diaferio for insightful discussions and
valuable comments on the draft. AL has been supported by the EU H2020-MSCA-
ITN-2019 Project 860744 ‘BiD4BESt: Big Data applications for Black Hole
Evolution Studies’ and by the PRIN MIUR 2017 prot. 20173ML3WW, ‘Opening the
ALMA window on the cosmic evolution of gas, stars, and supermassive black
holes’. BSH is supported by the INFN INDARK grant and acknowledges support
from the COSMOS project of the Italian Space Agency (cosmosnet.it). CB
acknowledges support from the COSMOS project of the Italian Space Agency
(cosmosnet.it), and the INDARK Initiative of the INFN
(web.infn.it/CSN4/IS/Linea5/InDark).

## References

  * Eckert  _et al._ (2022) D. Eckert, S. Ettori, E. Pointecouteau, R. F. J. van der Burg, and S. I. Loubser, Astronomy & Astrophysics 662, A123 (2022).
  * Aarseth  _et al._ (1979) S. J. Aarseth, I. Gott, J. R., and E. L. Turner, Astrophys. J. 228, 664 (1979).
  * de Andres  _et al._ (2022) D. de Andres, G. Yepes, F. Sembolini, G. Martí nez-Muñoz, W. Cui, F. Robledo, C.-H. Chuang, and E. Rasia, Monthly Notices of the Royal Astronomical Society 518, 111 (2022).
  * Cui  _et al._ (2018) W. Cui, A. Knebe, G. Yepes, F. Pearce, C. Power, R. Dave, A. Arth, and et al., Monthly Notices of the Royal Astronomical Society 480, 2898 (2018).
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Adam  _et al._ (2016) R. Adam _et al._ (Planck), Astron. Astrophys. 594, A1 (2016), arXiv:1502.01582 [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] .
  * Eckert  _et al._ (2016) D. Eckert _et al._ , Astron. Astrophys. 592, A12 (2016), arXiv:1512.03814 [astro-ph.CO] .
  * 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._ (2019a) 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 (2019a), arXiv:1805.00035 [astro-ph.CO] .
  * 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] .
  * Ettori  _et al._ (2019b) S. Ettori, V. Ghirardini, D. Eckert, E. Pointecouteau, F. Gastaldello, M. Sereno, M. Gaspari, S. Ghizzardi, M. Roncarelli, and M. Rossetti, Astronomy and Astrophysics 621, A39 (2019b).
  * Murray  _et al._ (2022) C. Murray, J. G. Bartlett, E. Artis, and J.-B. Melin, Monthly Notices of the Royal Astronomical Society 512, 4785 (2022).
  * Diaferio (1999) A. Diaferio, Monthly Notices of the Royal Astronomical Society 309, 610 (1999), https://academic.oup.com/mnras/article-pdf/309/3/610/2965424/309-3-610.pdf .
  * Biviano  _et al._ (2023) A. Biviano, L. Pizzuti, A. Mercurio, B. Sartoris, P. Rosati, S. Ettori, M. Girardi, C. Grillo, G. B. Caminha, and M. Nonino, The Astrophysical Journal 958, 148 (2023).
  * von der Linden  _et al._ (2014) A. von der Linden, M. T. Allen, D. E. Applegate, P. L. Kelly, S. W. Allen, H. Ebeling, P. R. Burchat, D. L. Burke, D. Donovan, R. G. Morris, R. Blandford, T. Erben, and A. Mantz, Monthly Notices of the Royal Astronomical Society 439, 2 (2014).
  * Hurier and Angulo (2018) G. Hurier and R. E. Angulo, Astronomy & Astrophysics 610, L4 (2018).
  * Maughan  _et al._ (2016) B. J. Maughan, P. A. Giles, K. J. Rines, A. Diaferio, M. J. Geller, N. V. D. Pyl, and M. Bonamente, Monthly Notices of the Royal Astronomical Society 461, 4182 (2016).
  * Umetsu (2020) K. Umetsu, Astron. Astrophys. Rev. 28, 7 (2020), arXiv:2007.00506 [astro-ph.CO] .
  * et al (2018) A. A. et al, Monthly Notices of the Royal Astronomical Society 477, 4285 (2018).
  * et al. (2018) O. T. et al., Monthly Notices of the Royal Astronomical Society 477, 5657 (2018).
  * Gifford  _et al._ (2013) D. Gifford, C. Miller, and N. Kern, Astrophys. J. 773, 116 (2013), arXiv:1307.0017 [astro-ph.CO] .
  * Serra  _et al._ (2010) A. L. Serra, A. Diaferio, G. Murante, and S. Borgani, Monthly Notices of the Royal Astronomical Society , no (2010).
  * et al. (2023) Z. P. et al., “A measurement of gravitational lensing of the cosmic microwave background using spt-3g 2018 data,” (2023), arXiv:2308.11608 [astro-ph.CO] .
  * et al. (2012) F. W. H. et al., The Astrophysical Journal 758, 68 (2012).
  * Diaferio  _et al._ (2005) A. Diaferio, M. J. Geller, and K. J. Rines, The Astrophysical Journal 628, L97 (2005).
  * Diaferio and Geller (1997) A. Diaferio and M. J. Geller, The Astrophysical Journal 481, 633 (1997).
  * Diaferio (2009) A. Diaferio, arXiv e-prints , arXiv:0901.0868 (2009), arXiv:0901.0868 [astro-ph.CO] .
  * Logan  _et al._ (2022) C. H. A. Logan, B. J. Maughan, A. Diaferio, R. T. Duffy, M. J. Geller, K. Rines, and J. Sohn, Astronomy & Astrophysics 665, A124 (2022).
  * Andreon  _et al._ (2017) S. Andreon, G. Trinchieri, A. Moretti, and J. Wang, Astronomy &Astrophysics 606, A25 (2017).
  * 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] .
  * Babichev and Deffayet (2013) E. Babichev and C. Deffayet, Classical and Quantum Gravity 30, 184001 (2013).
  * Khoury and Weltman (2004a) J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004a), arXiv:astro-ph/0309300 .
  * Zaregonbadi  _et al._ (2022) R. Zaregonbadi, N. Saba, and M. Farhoudi, Eur. Phys. J. C 82, 730 (2022), arXiv:2207.14475 [gr-qc] .
  * Ivanov and Wellenzohn (2020) A. N. Ivanov and M. Wellenzohn, Universe 6, 221 (2020), arXiv:1607.00884 [gr-qc] .
  * Kraiselburd  _et al._ (2018) L. Kraiselburd, S. J. Landau, M. Salgado, D. Sudarsky, and H. Vucetich, Phys. Rev. D 97, 104044 (2018), arXiv:1511.06307 [gr-qc] .
  * Tsujikawa  _et al._ (2009) S. Tsujikawa, T. Tamaki, and R. Tavakol, JCAP 05, 020 (2009), arXiv:0901.3226 [gr-qc] .
  * Vainshtein (1972) A. I. Vainshtein, Phys. Lett. B 39, 393 (1972).
  * 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, Physical Review D 97 (2018), 10.1103/physrevd.97.021301.
  * Langlois  _et al._ (2018) D. Langlois, R. Saito, D. Yamauchi, and K. Noui, Physical Review D 97 (2018), 10.1103/physrevd.97.061501.
  * Bartolo  _et al._ (2018) N. Bartolo, P. Karmakar, S. Matarrese, and M. Scomparin, Journal of Cosmology and Astroparticle Physics 2018, 048 (2018).
  * Dima and Vernizzi (2018) A. Dima and F. Vernizzi, Physical Review D 97 (2018), 10.1103/physrevd.97.101302.
  * Hirano  _et al._ (2019) S. Hirano, T. Kobayashi, D. Yamauchi, and S. Yokoyama, (2019), arXiv:1902.02946 [astro-ph.CO] .
  * Crisostomi  _et al._ (2019) M. Crisostomi, M. Lewandowski, and F. Vernizzi, Phys. Rev. D 100, 024025 (2019), arXiv:1903.11591 [gr-qc] .
  * Terukina  _et al._ (2014) A. Terukina, L. Lombriser, K. Yamamoto, D. Bacon, K. Koyama, and R. C. Nichol, J. Cosmology Astropart. Phys. 2014, 013 (2014).
  * Haridasu  _et al._ (2023) B. S. Haridasu, P. Karmakar, M. De Petris, V. F. Cardone, and R. Maoli, Phys. Rev. D 107, 124059 (2023), arXiv:2111.01101 [astro-ph.CO] .
  * Kaiser (1987) N. Kaiser, Mon. Not. Roy. Astron. Soc. 227, 1 (1987).
  * van Haarlem and van de Weygaert (1993) M. van Haarlem and R. van de Weygaert, Astrophys. J. 418, 544 (1993).
  * Regos and Geller (1989) E. Regos and M. J. Geller, AJ 98, 755 (1989).
  * Mamon  _et al._ (2013) G. A. Mamon, A. Biviano, and G. Boué, Mon. Not. Roy. Astron. Soc. 429, 3079 (2013), arXiv:1212.1455 [astro-ph.CO] .
  * Read  _et al._ (2021) J. I. Read, G. A. Mamon, E. Vasiliev, L. L. Watkins, M. G. Walker, J. Peñarrubia, M. Wilkinson, W. Dehnen, and P. Das, Mon. Not. Roy. Astron. Soc. 501, 978 (2021), arXiv:2011.09493 [astro-ph.GA] .
  * Binney and Mamon (1982) J. Binney and G. A. Mamon, Mon. Not. Roy. Astron. Soc. 200, 361 (1982).
  * Mamon  _et al._ (2019) G. A. Mamon, A. Cava, A. Biviano, A. Moretti, B. Poggianti, and D. Bettoni, Astron. Astrophys. 631, A131 (2019), arXiv:1901.06393 [astro-ph.GA] .
  * Geller  _et al._ (1999) M. J. Geller, A. Diaferio, and M. J. Kurtz, Astrophys. J. Lett. 517, L23 (1999), arXiv:astro-ph/9903305 .
  * Biviano and Girardi (2003) A. Biviano and M. Girardi, The Astrophysical Journal 585, 205 (2003).
  * Khoury and Weltman (2004b) J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004b).
  * Khoury (2013) J. Khoury, “Les houches lectures on physics beyond the standard model of cosmology,” (2013), arXiv:1312.2006 [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] .
  * Paillas  _et al._ (2019) E. Paillas, M. Cautun, B. Li, Y.-C. Cai, N. Padilla, J. Armijo, and S. Bose, Monthly Notices of the Royal Astronomical Society 484, 1149 (2019).
  * Brando  _et al._ (2023) G. Brando, K. Koyama, and H. A. Winther, Journal of Cosmology and Astroparticle Physics 2023, 045 (2023).
  * Cardone  _et al._ (2020) V. F. Cardone _et al._ , arXiv (2020), 10.1103/PhysRevD.103.064065, 2008.00580 .
  * Langlois  _et al._ (2017) D. Langlois, M. Mancarella, K. Noui, and F. Vernizzi, Journal of Cosmology and Astroparticle Physics 2017, 033 (2017).
  * 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] .
  * Haridasu  _et al._ (2021) B. S. Haridasu, P. Karmakar, M. De Petris, V. F. Cardone, and R. Maoli, in  _mm Universe @ NIKA2_ (2021) arXiv:2111.01720 [astro-ph.CO] .
  * Gandolfi  _et al._ (2023) G. Gandolfi, B. S. Haridasu, S. Liberati, and A. Lapi, (2023), arXiv:2305.13974 [astro-ph.CO] .
  * Benetti  _et al._ (2023) F. Benetti, A. Lapi, G. Gandolfi, B. S. Haridasu, and L. Danese, Universe 9, 329 (2023), arXiv:2307.04655 [astro-ph.CO] .
  * Sohn  _et al._ (2017) J. Sohn, M. J. Geller, H. J. Zahid, D. G. Fabricant, A. Diaferio, and K. J. Rines, The Astrophysical Journal Supplement Series 229, 20 (2017).
  * Liu  _et al._ (2018) A. Liu, H. Yu, A. Diaferio, P. Tozzi, H. S. Hwang, K. Umetsu, N. Okabe, and L.-L. Yang, The Astrophysical Journal 863, 102 (2018).
  * Boumechta  _et al._ (2023) Y. Boumechta, B. S. Haridasu, L. Pizzuti, M. A. Butt, C. Baccigalupi, and A. Lapi, Phys. Rev. D 108, 044007 (2023), arXiv:2303.02074 [astro-ph.CO] .
  * Schaller  _et al._ (2015) M. Schaller, C. S. Frenk, R. G. Bower, T. Theuns, J. Trayford, R. A. Crain, M. Furlong, J. Schaye, C. Dalla Vecchia, and I. G. McCarthy, Monthly Notices of the Royal Astronomical Society 452, 343 (2015), https://academic.oup.com/mnras/article-pdf/452/1/343/4932119/stv1341.pdf .
  * Hogan  _et al._ (2017) M. T. Hogan, B. R. McNamara, F. Pulido, P. E. J. Nulsen, H. R. Russell, A. N. Vantyghem, A. C. Edge, and R. A. Main, The Astrophysical Journal 837, 51 (2017).
  * Sartoris  _et al._ (2020) B. Sartoris, A. Biviano, P. Rosati, A. Mercurio, C. Grillo, S. Ettori, M. Nonino, K. Umetsu, P. Bergamini, G. B. Caminha, and M. Girardi, Astron. Astrophys. 637, A34 (2020), arXiv:2003.08475 [astro-ph.CO] .
  * 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] .
  * Wing and Blanton (2011) J. D. Wing and E. L. Blanton, The Astronomical Journal 141, 88 (2011).
  * Pizzardo  _et al._ (2023a) M. Pizzardo, M. J. Geller, S. J. Kenyon, I. Damjanov, and A. Diaferio, Astron. Astrophys. 675, A56 (2023a), arXiv:2303.18147 [astro-ph.CO] .
  * Springel  _et al._ (2018) V. Springel _et al._ , Mon. Not. Roy. Astron. Soc. 475, 676 (2018), arXiv:1707.03397 [astro-ph.GA] .
  * Herbonnet  _et al._ (2020) R. Herbonnet, C. Sifon, H. Hoekstra, Y. Bahe, R. F. J. van der Burg, J.-B. Melin, A. von der Linden, D. Sand, S. Kay, and D. Barnes, Mon. Not. Roy. Astron. Soc. 497, 4684 (2020), arXiv:1912.04414 [astro-ph.CO] .
  * Sakstein (2015) J. Sakstein, Phys. Rev. D92, 124045 (2015), arXiv:1511.01685 [astro-ph.CO] .
  * Pizzuti  _et al._ (2020) L. Pizzuti, I. D. Saltas, and L. Amendola, (2020), arXiv:2011.15089 [astro-ph.CO] .
  * Pizzuti  _et al._ (2023) L. Pizzuti, I. D. Saltas, A. Biviano, G. Mamon, and L. Amendola, J. Open Source Softw. 8, 4800 (2023), arXiv:2201.07194 [astro-ph.CO] .
  * Gifford  _et al._ (2017) D. Gifford, N. Kern, and C. J. Miller, Astrophys. J. 834, 204 (2017), arXiv:1612.05562 [astro-ph.CO] .
  * Svensmark  _et al._ (2015) J. Svensmark, R. Wojtak, and S. H. Hansen, Mon. Not. Roy. Astron. Soc. 448, 1644 (2015), arXiv:1405.0284 [astro-ph.CO] .
  * Pizzardo  _et al._ (2023b) M. Pizzardo, M. J. Geller, S. J. Kenyon, I. Damjanov, and A. Diaferio, (2023b), arXiv:2307.07398 [astro-ph.CO] .
  * Wilcox  _et al._ (2016) H. Wilcox, R. C. Nichol, G.-B. Zhao, D. Bacon, K. Koyama, and A. K. Romer, Mon. Not. Roy. Astron. Soc. 462, 715 (2016), arXiv:1603.05911 [astro-ph.CO] .

##  Appendix A Error estimation in the caustic technique

In this section, we highlight a few systematics that affect the estimated
caustic mass. The caustic technique heavily relies on the methods utilized to
set the so-called ‘Caustic’ surface. The caustic surface is estimated as the
edge of the galaxy cluster assessed through the distribution of the galaxies
in the line-of-sight velocity and the projected radius from the center of the
clusters. We utilize our importance sampling-like approach instead of a
likelihood-based technique as the error estimation of the caustic surface
relies on the density of the galaxies, which can be estimated at every radius
instead of data with error bars that can be incorporated into a simple
Gaussian likelihood.

##  Appendix B Comment on the anisotropy profile $g_{\beta}(r)$ assumption

As mentioned in the main text, we have assumed a fixed anisotropy profile of
the galaxy phase space in the current analysis. To check the robustness of our
statement, we also fit the anisotropy profile using the MG-MAMPOSSt code
Pizzuti  _et al._ (2020), which provides joint mass density and orbit
anisotropy reconstructions in GR and several modified gravity scenarios. More
in detail, MG-MAMPOSSt relies on data in the projected phase space of member
galaxies; it solves the spherically-symmetric Jeans’ equation, assuming a
shape for the three-dimensional velocity distribution and dynamical relaxation
of the system. The code implements a parametric determination of the kinematic
quantities (i.e. gravitational potential, number density of galaxies, and
velocity anisotropy) by performing a Maximum Likelihood fit in the projected
phase space.

Assuming both GR and one of the modified gravity frameworks, namely Chameleon
screening, we reconstruct the mass distribution and the anisotropy profile of
A2029 and A2142. In all cases, we adopt an NFW model for the matter density
distribution, while the velocity anisotropy is parametrized with a generalized
Tiret model (e.g. Mamon  _et al._ (2019)),

$\beta_{\text{gT}}(r)=\beta_{0}+(\beta_{\infty}+\beta_{0})\frac{r}{r-r_{\beta}}\,$
(24)

where $\beta_{0}$, $\beta_{\infty}$ represent the central anisotropy and the
anisotropy at a large distance from the center, respectively; $r_{\beta}$ is a
characteristic radius which we assume to be equivalent to the scale radius
$r_{\text{s}}$ of the total mass profile. Equation 24 provides a quite general
description of possible orbit anisotropy profiles in galaxy clusters. For both
clusters, we consider galaxies that lie in projection within $R\sim
1.1\,R_{200}^{\text{HS}}$, to ensure the validity of Jean’s equation.
Nevertheless, we checked that reasonable changes ($10\%$) in that limit
provide negligible effects on the final results. We further exclude the
central region $R<0.05\,{\rm Mpc})$, where the dynamics are dominated by the
Brighter Central Galaxy.

We use than the anisotropy profile reconstructed with MG-MAMPOSSt to validate
the statement that our assumption has a negligible effect on the inferred
final caustic mass. Essentially, this is because the anisotropy profile, which
requires the knowledge of $R_{200}$, is degenerate with the cluster viral
radius, and hence with the total mass. However, since the hydrostatic
technique constraints the same quantity much more tightly than the anisotropy
profile itself, we find that fixing it to a constant is largely valid and
makes minimal changes to the conclusions made here.

##  Appendix C Incorrect assumptions of $\mathcal{F}_{\beta}(r)$ in modified
gravity

It is well established that the $\mathcal{F}_{\beta}(r)$ in the inner regions
of the cluster is not a constant and assuming so can lead to biases. In this
section, we show as an example, the potential of the caustic technique to very
well constrain the modified gravity scenarios in case the
$\mathcal{F}_{\beta}(r)$ derived from hydrostatic equilibrium differs from
that constrained using alternate techniques or simulations. This is
particularly important in the modified gravity scenarios, where the
constraints on the mass are relaxed when using the hydrostatic equilibrium
data (see Figure 9). As a test case, we show this through the example of the
cluster A2029, when assuming $\mathcal{F}_{\beta}(r)=0.5$ as an independent
estimate. When performing the joint analysis the essential improvement in the
constraints on the modified gravity model parameters is valid if and only if
the caustic masses are inferred utilizing the apriori constrained masses using
the hydrostatic technique. For instance, in Figure 8, we show the joint
constraints from the importance sampling against the hydrostatic-only
constraints, wherein the caustic technique assumes a constant
$\mathcal{F}_{\beta}(r)=0.5$. One can immediately notice the well-bounded
constraints on the parameters $\\{\phi_{\infty,2},\beta_{2}\\}$, which are
spurious posteriors indicating evidence for modified gravity, in the current
comparison. However, as shown in the Figures 6 and 7, for the cluster A2029
while the caustic technique helps to tighten the posterior assuming the
correct $\mathcal{F}_{\beta}(r)$, does not indicate any deviation from GR.
Needless to say, we make equivalent conclusions also in the case of the
cluster A2142.

Figure 8: Left: We show the constraints on the modified gravity parameters in
the case of the chameleon screening when assuming a
$\mathcal{F}_{\beta}(r)=0.5$. Right: Same as Left panel, but for the case of
the Vainshtein screening.

##  Appendix D Mass profiles

In Figure 9, we show the mass profiles in various scenarios. We find that the
model-independent caustic mass profile is in very good agreement with the NFW-
based mass profile constrained using the hydrostatic data. We also find an
overall slightly lower mass profile when assuming the
$\mathcal{F}_{\beta}(r)=0.5$ and a very good agreement when using
$\mathcal{F}_{\beta}(r)=0.65$, while showing no discernable change in the
shape of the profile. The caustic mass profile also shows a tendency to get
flattened in the outskirts of the cluster, which is a better assessment of the
mass profile in contrast to assuming a mass profile such as NFW which
asymptotically increases. While we retain the comparison in the main text to
the mass $M_{500}$ of the cluster, which is convenient for performing
statistical analysis, we could also place upper limits on the total mass of
the cluster using the caustic mass profile. Note that the mass profiles
reconstructed using the caustic technique flatten out, seemingly suggesting an
asymptotically converging total mass of the cluster. However, this is only an
artifact of the reconstruction technique, wherein the caustic surface,
$\mathcal{A}^{2}(r)$ artificially goes to zero where no galaxies are present.
This in turn implies the integrand in Equation 3 goes to zero providing
constant mass values after a certain radial range. Therefore, we note that
while we show the mass profiles going out to 6 Mpc, the nature of the increase
of mass should be inferred from the NFW-based extrapolation or the flattened
caustic mass profiles. In this context, it is also possible to utilize a value
of mass say at $\Delta=100$ to perform our joint analysis, which we intend to
leave for a future investigation while improving the joint analysis.

Figure 9: We compare the mean and the $68\%$ confidence level regions of the
mass profiles obtained using the hydrostatics and caustics techniques. Top:
Mass profiles in the case of GR (left), chameleon screening (centre) and
Vainshtein screening (right). In each of the panels, we show the NFW-based
mass profile (yellow), model-independent reconstructed mass profile using the
caustic technique, assuming hydrostatic $\mathcal{F}_{\beta}(r)$ (green) and
constant $\mathcal{F}_{\beta}(r)=0.65$ (red). Bottom: Same as Top panel, for
the cluster A2142.