are required than for standard CDM runs.

For the ${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ SIDM merging run it
was found after some testing that the best results are obtained for
$r_{c}^{\rm{prm}}\sim 400\rm{\,kpc}$ and cluster gas fractions
$(f_{g1},f_{g2})=(0.12,0.12)$. We label this model as XDBf_sa. Similarly, for
the ${\sigma_{DM}}/{m_{X}}=4\rm{\,cm^{2}}{gr^{-1}}$ simulation it was found
that $r_{c}^{\rm{prm}}\sim 300\rm{\,kpc}$ and $(f_{g1},f_{g2})=(0.12,0.14)$
are the optimal settings to match the observed X-ray morphology. We label this
model as XDBf_sb. We report in Table 7 the most relevant initial merger
parameters of the two models; in both cases the initial gas density profile of
the secondary is the same as in the standard CDM runs. Figure 12 shows the
mock X-ray maps extracted from these two SIDM merger simulations.

For the two merging clusters, a specific signature of SIDM is the expected
offsets of their DM centroids from those of the other cluster mass components.
For this reason, in each panel of Figure 12 we report for each simulation the
relative distance between the different centroids.

Specifically, for the SE cluster we indicate the distance of the BCG to the
the DM centroid ($d_{\rm{BCG-DM}}$), the BCG to the X-ray peak
($d_{\rm{BCG-X}}$) and the DM to the X-ray peak ($d_{\rm{X-DM}}$). For the NW
cluster we report only the distance of the SZ to the DM centroid
($d^{NW}_{\rm{SZ-DM}}$), but not that of the DM to the X-ray peak because the
gas structure of the primary has been destroyed during the collision.

In order to assess how well the SIDM simulations can reproduce the observed
offsets, we also show in Figure 12 the measured positions of the different
centroids. Data points are extracted from the two panels of Figure 6 of Kim et
al. (2021), the top (bottom) panel is for the NW (SE) cluster. For the sake of
comparison we report here for the two clusters the relative location of
several peaks, as measured with respect the positions
$(\Delta_{x},\Delta_{y})=(0.0^{\prime\prime},0.0^{\prime\prime})$ of the
cluster mass centroids. For the SE cluster the measured relative positions are
: ${\bf\Delta}^{\rm{SE-BCG}}=(-6.8^{\prime\prime},0^{\prime\prime})$ and
${\bf\Delta}^{\rm{SE-X}}=(-9.6^{\prime\prime},-8.1^{\prime\prime})$; while for
the NW cluster: ${\bf\Delta}^{\rm{NW-
SZ}}=(-16^{\prime\prime},-16.7^{\prime\prime})$ and ${\bf\Delta}^{\rm{NW-
BCG}}=(-7.15^{\prime\prime},11.4^{\prime\prime})$.

The NW cluster does not possess a BCG, but its galaxy number density peak
shows a large offset with respect the mass centroid position (Jee et al.
2014). With the notation ${\bf\Delta}^{\rm{NW-BCG}}$ we then indicate the
relative position of the galaxy number peak as shown in the top panel of
Figure 6 of Kim et al. (2021). The measured relative offsets are indicated in
Figure 12 with filled circles, their positions are with respect that of the DM
centroids and in accord with the adopted cosmology we set the spatial
conversion of angular coordinates to $1^{\prime}\simeq 470\rm{\,kpc}$.

The absence of observational errors for the measured offsets is clearly a
critical issue in order to set constraints on DM cross-sections as well as to
rule out the null hypothesis of zero offsets. According to Kim et al. (2021),
the statistical significance of the X-ray peak offset from the SE mass
centroid is at $\sim 2\sigma$ level. Similarly, the galaxy number density peak
is offset from the NW (SE) mass centroid at $\sim 2\sigma$ ($\sim 6.6\sigma$)
level. The statistical significance of the SZ peak offset from the NW mass
peak is at $\sim 3.7\sigma$ level, although Kim et al. (2021) caution that the
interpretation of this offset is complicated by the low angular resolution of
the beam size ($\sim 1^{\prime}.4$). Overall, these thresholds allow us to
exclude with high significance zero size offsets.

For the two SIDM merging simulations the X-ray maps of Figure 12 show the
spatial locations of the different centroids, as extracted from the
simulations. Their relative positions constitute one of the most interesting
results of our study, and we suggest that these findings pose a severe
challenge to the collisionless standard CDM model.

As already outlined, the most striking feature is the position for the SE
cluster of the X-ray peak, which is now leading the DM mass centroid.
Moreover, the relative distance $d_{\rm{X-DM}}$ between the two centroids
varies from $d_{\rm{X-DM}}\sim 60\rm{\,kpc}$
(${\sigma_{DM}}/{m_{X}}=4\rm{\,cm^{2}}{gr^{-1}}$) up to $d_{\rm{X-DM}}\sim
120\rm{\,kpc}$ (${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$). The
agreement of this range of values with the measured offset $d^{SE}_{X-DM}\sim
100\rm{\,kpc}$ can be considered significant, given the observational
uncertainties (see later).

The offset of the SZ peak is similarly affected in an SIDM merging scenario.
For the standard CDM merging simulations of Figure 10, the distance of the SZ
peak from the NW DM centroid is about $d^{NW}_{SZ-DM}\sim 330\rm{\,kpc}$. This
value is similar to that found by ZYL15 for their model B and higher than the
measured value of about $\sim 150\rm{\,kpc}$ (Jee et al. 2014). However for
the SIDM simulations of Figure 12 this distance is now reduced to $d^{NW}_{SZ-
DM}\sim 230\rm{\,kpc}$, in better agreement with the data.

Another important result which follows from an SIDM scenario in major cluster
mergers is the expected presence of BCG-DM offsets, with BCG centroids that
exhibit miscentred positions with respect to the centers of their parent DM
halos. Figure 12 shows that according to value of ${\sigma_{DM}}/{m_{X}}$, the
BCG-DM offsets for the SE cluster lie in the range $d_{\rm{BCG-DM}}\sim
60-120\rm{\,kpc}$, which agrees fairly well with the observed offset of about
$\sim 60\rm{\,kpc}$. Note that for the NW cluster, the position of the
simulated BCG strictly coincides with that of the observed peak galaxy number
density.

In SIDM mergers, positive BCG-DM offsets are expected after the first
pericenter passage. This effect can be interpreted as a result of the
deceleration of the DM component, this in turn is due to the exchange of
energy between the DM halos which occurs during the cluster collision. As
mentioned in the discussion of the results from the
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ SIDM simulations of Figure 11,
positive BCG-DM offsets are a transient phenomena. As a BCG starts to climb
out of the cluster potential after the pericenter passage, it will experience
a gravitational pull (Harvey et al. 2014) from the DM mass component which
will reduce its forward momentum as well as the size of the offset.

This restoring force has the effect of reducing the BCG bulk velocities in the
SIDM mergers of Figure 12, similarly reducing the relative mean radial
velocity $V^{s}_{r}$ along the line-of-sight between the two SE and NW BCG
components. As a result the BCG star velocities can no longer be considered a
fair proxy for the DM particle velocities. This bias strongly helps to reduce
the discrepancy between the mean relative velocities $V_{r}$, extracted from
the DM particle velocities of the fiducial models of Section 3.1 (see Figure
7), and those estimated from galaxy-based spectroscopic measurements
(Menanteau et al. 2012).

In Figure 12 the measured velocity offsets $V_{r}^{s}$ between the two BCGs
are reported for the two merger models. The values of $V_{r}^{s}$ range from
$V_{r}^{s}\sim 300\rm{\,km\leavevmode\nobreak\ s^{-1}}$ for
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ up to $V_{r}^{s}\sim
650\rm{\,km\leavevmode\nobreak\ s^{-1}}$ when
${\sigma_{DM}}/{m_{X}}=4\rm{\,cm^{2}}{gr^{-1}}$. These values can be
contrasted with that given by Menanteau et al. (2012), who reported
$V_{r}^{s}=598\pm 96\rm{\,km\leavevmode\nobreak\ s^{-1}}$ for the relative
radial velocity of the SE cluster with respect the NW component. We therefore
conclude that the SIDM merger models of the El Gordo cluster presented here,
predict mean relative radial velocities between their stellar components in
much better agreement with data than in the collisionless CDM cases.

Although the simulations of the two SIDM merger models are undoubtedly able to
reproduce the observed offsets and the relative cluster peculiar velocity,
nonetheless the two X-ray maps of Figure 12 exhibit a twin-tailed X-ray
morphology much less defined than that displayed by the merger model Bf_rc20
in Figure 2. In fact, for model XDBf_sa the tails behind the secondary are
almost absent.

Moreover, the X-ray emission in the outer regions behind the secondary is
significantly less than in model Bf_rc20. This can be clearly seen by the the
location of the third contour level ($2.9\cdot
10^{-1}\rm{\,counts\leavevmode\nobreak\ arcsec^{-2}}$ from the inside), which
in the X-ray map of model Bf_rc20 ( top right panel of Figure 2) has a minimum
along the line joining the two DM centroids at about $\sim 600\rm{\,kpc}$ from
the X-ray peak of the secondary. For the merger simulation of model XDBf_sb
this minimum is much less pronounced and at about $\sim 400\rm{\,kpc}$ from
the centroid of the X-ray emission.

These shortcomings of the two SIDM merger models in reproducing the observed
X-ray morphology are present even after the adoption of initially higher gas
fractions and of larger gas scale radii for the primary. These initial
condition parameters were purposely introduced to compensate for the much
weaker resiliency of the gas structures seen in the SIDM merging runs of
Figure 11.

We argue that this difficulty of SIDM merging simulations in matching the
observed X-ray morphology is strictly related to the complex dynamical
interplay which occurs during the cluster collision between baryons and SIDM.
We suggest that in order to extract from SIDM simulations X-ray maps more in
agreement with observations, it is necessary to explore a much wider parameter
space for the initial gas density profiles of the two clusters. The results
presented here for these simulations should then be regarded as preliminary.
We postpone to a future work a much more thorough study on the impact of SIDM
on the X-ray morphology in cluster mergers.

##  4 Summary and Conclusions

We have constructed a large set of N-body/hydrodynamical binary cluster merger
simulations, aimed at studying the merging cluster El Gordo. Specifically, the
main goal of this paper is to assess whether it is possible for the merger
simulations to reproduce the observed twin-tailed X-ray morphology, while
satisfying other observational constraints.

To this end we considered a wide range of initial conditions, aimed at
performing both off-axis and head-on merging runs. The initial conditions of
the binary merging clusters are constrained to be consistent with recent
lensing-based mass measurements (Diego et al. 2020; Kim et al. 2021), so that
we impose for the mass of the primary the following range of values:
$10^{15}\rm{\,M_{\odot}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}M\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1.6\cdot 10^{15}\rm{\,M_{\odot}}$.
A summary of our main results is as follows:

(i) The observed twin-tailed X-ray morphology, as well as other observational
constraints, are well matched by the off-center fiducial merger models
Bf_rc20, Bg_rc20 and Bl_rc24 of Section 3.1 that have collision velocities and
impact parameters in the range $2,000\rm{\,km\leavevmode\nobreak\
s^{-1}}\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}V\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2,500\rm{\,km\leavevmode\nobreak\
s^{-1}}$ and $600\rm{\,kpc}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}P\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}800\rm{\,kpc}$, respectively.

These findings demonstrate that in simulations with cluster masses lower than
previously considered (Donnert 2014; Molnar & Broadhurst 2015; Zhang et al.
2015, 2018), it is possible to consistently reproduce the observed twin-tailed
X-ray morphology seen in the cluster El Gordo. Finally, the observational
constraints discussed in Section 3.1 tend to favor merger models with a mass
of the primary of about $M\simeq 1.6\cdot 10^{15}\rm{\,M_{\odot}}$.

(ii) One of the most interesting aspects of the galaxy cluster El Gordo is the
odd location in the SE cluster of the X-ray peak, which in the outgoing
scenario seems to be ahead of the DM lensing peak. To resolve this issue a
returning scenario has been invoked (Ng et al. 2015), so that the merging
cluster is then observed in a post-apocenter phase and the X-ray peak would be
trailing the DM peak.

This scenario has been extensively investigated in Section 3.2, in which we
constructed a large set of head-on merging collisions aimed at reproducing the
observed twin-tailed X-ray morphology of the El Gordo cluster. For the adopted
initial collision parameter space the results of Figure 9 demonstrate that a
returning scenario should be considered untenable, as the morphology of all of
the X-ray maps extracted from the simulations appears inconsistent with the
observed X-ray morphology. We argued in Section 3.2 that this behavior is
strictly related to the lifetime of the post-pericenter X-ray structures,
which turned out to be much less than the elapsed time needed for the DM
component of the SE cluster to reach the apocenter and return.

We conclude that the likelihood of head-on collisions being able to match the
observational constraints for the merging cluster El Gordo is very low, given
the wide range of initial merging parameter space that has been considered.

(iii) We resimulated two of the fiducial models (Bf_rc20 and Bl_rc24) of
Section 3.1 by adding to the initial mass components of each of the two halos
a star particle distribution, with the purpose of representing the BCG’s mass
contribution to the initial halo profiles. The resulting images for these two
runs (models DBf_rc20 and DBl_rc24) are shown in Figure 10. The positions of
the centroids of the different mass components clearly indicate that, in both
of the mergers, the collision between the two clusters has not produced any
appreciable offset between the positions of the BCG centroids relative to that
of the DM halos.

In a collisionless CDM scenario this leaves us with the problem of explaining
the presence and the magnitude of the BCG to DM offset seen in the El Gordo
cluster. It must be emphasized that galaxy-DM and BCG-DM offsets are relative
rare but have been measured in other merging clusters, (see, for example,
Table 2 of Kim et al. 2017).

According to Martel et al. (2014) such offsets are produced as a result of
violent cluster collisions which strongly perturb the clusters and leave them
in a non-equilibrium state, with the BCGs being displaced from the cluster
centers. While in principle this scenario cannot be ruled out for explaining
the observed BCG to DM offset seen the for El Gordo SE cluster, nonetheless we
argue that it appears difficult for the gas structure of the SE cluster to
maintain its integrity after a cluster collision which must be at the same
time violent enough to perturb the BCG from its original position.

Finally, it should be emphasized that the observed BCG to DM offset
constitutes another serious shortcoming for the returning scenario: as shown
by the the simulations of Figure 10, this offset would remain unexplained in
such a scenario.

(iv) Motivated by these findings we studied in Section 3.3 results extracted
from SIDM simulations of the El Gordo cluster, in which DM interactions are
described by a simple model with a velocity independent elastic DM cross-
section.

The most important striking feature that emerges from the SIDM merging models
of the El Gordo cluster is that the observed offsets between the different
mass components are well reproduced by the simulations XDBf_sa and XDBf_sb of
Table 7. These simulations have a DM cross-section in the range
${\sigma_{DM}}/{m_{X}}\sim 4-5\rm{\,cm^{2}}{gr^{-1}}$ and initial collision
parameters of model Bf_rc20 in Table 2.

Another significant aspect of these merger models is the value of the mean
relative line-of-sight radial velocity between the two BCGs, which in contrast
with previous runs (see Figure 7) is now of the order of several hundreds of
$\rm{\,km\leavevmode\nobreak\ s^{-1}}$ and in much better agreement with data.
As discussed in Section 3.3, this is a direct consequence of the gravitational
pull experienced by the two BCGs because of the dissipative interactions now
present between the two DM halos.

Although we found it difficult for these models to reproduce the observed
X-ray morphology as well as in the standard CDM runs of Section 3.1,
nonetheless we argue that these findings provide significant support to the
possibility that in the El Gordo cluster the DM behavior exhibits collisional
properties.

Additionally, it must be stressed that a significant advantage of the SIDM
merging models is to offer a natural explanation for all of the observed
offsets, whilst in the standard CDM scenario the observed separations between
the different centroids remain either unexplained (X-ray peak to DM) or
fraught with difficulties (BCG to DM).

From these conclusions follow two main consequences that need to be discussed.
The first is that the two SIDM models of Table 7 do not admit free parameters;
the observed offsets are well reproduced only for a DM scattering cross-
section in the range ${\sigma_{DM}}/{m_{X}}\sim 4-5\rm{\,cm^{2}}{gr^{-1}}$ and
initial collision parameters given by model Bf_rc20. This is a direct result
of the post-pericenter evolution in an SIDM scenario, for which the observed
separations between the different centroids are locked to the measured
distance $d_{DM}$ between the DM mass centroids, which must be about
$d_{DM}\sim 700\rm{\,kpc}$ at the present epoch.

As a result, these observational constraints cannot be satisfied by simply
rescaling the value of ${\sigma_{DM}}/{m_{X}}$, if one assumes lower masses
for the merging cluster than in model Bf_rc20:
$(M^{(1)}_{200},\leavevmode\nobreak\ M^{(2)}_{200})=(1.6,0.69)\cdot
10^{15}\rm{\,M_{\odot}}$. This value of the primary’s mass is higher than that
obtained from recent lensing-based estimates (see Table 2 of Kim et al. 2021),
but still within the confidence limits. If more accurate measurements were to
yield tighter limits, clearly the SIDM merger models presented here could be
in jeopardy.

The other striking feature of the SIDM merger models XDBf_sa and XDBf_sb is
that the best match to the data is obtained only for values of
${\sigma_{DM}}/{m_{X}}$ in the very narrow range between $\sim
4\rm{\,cm^{2}}{gr^{-1}}$ and $\sim 5\rm{\,cm^{2}}{gr^{-1}}$. Such values are
strongly excluded by the present upper limits on the self-interaction DM
cross-section, which on galaxy cluster scales have been derived independently
by various authors (Robertson et al. 2017, 2019; Kim et al. 2017; Wittman et
al. 2018; Harvey et al. 2019; Shen et al. 2022; Cross et al. 2023). These
upper limits are at best of order unity (${\sigma_{DM}}/{m_{X}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1\rm{\,cm^{2}}{gr^{-1}}$), although
Wittman et al. (2018) claim that the constraint can be relaxed up to
${\sigma_{DM}}/{m_{X}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2\rm{\,cm^{2}}{gr^{-1}}$ by
properly averaging galaxy-DM offsets in a sample of merging clusters.

For example in the case of the Bullet Cluster the absence, within the
observational errors on the centroid positions, of an observed galaxy-DM
separation ($\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}20\rm{\,kpc}$,
Randall et al. 2008) has been used in SIDM merging simulations (Markevitch et
al. 2004; Randall et al. 2008; Robertson et al. 2017) to constrain the allowed
values of the DM cross-section to lie below ${\sigma_{DM}}/{m_{X}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1\rm{\,cm^{2}}{gr^{-1}}$. This
evidently implies that an SIDM merging simulation of the Bullet Cluster, with
the values derived here for the DM cross-section (${\sigma_{DM}}/{m_{X}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}4-5\rm{\,cm^{2}}{gr^{-1}}$), would
lead to simulation results clearly incompatible with the observational
constraints.

However, it must be emphasized the difficulty to quantitatively assess the
statistical significance of the above findings. This is because of the absence
of observational errors for the various offsets derived from Figure 6 of Kim
et al. (2021).

Kim et al. (2017) presented a set of SIDM merging simulations of equal-mass
massive clusters, with masses of the order of $\sim 10^{15}\rm{\,M_{\odot}}$.
They contrasted the expected galaxy-DM offsets from the SIDM simulations
against measured galaxy-DM offsets in observed mergers. The authors were able
to derive an upper limit on the galaxy-DM offset of $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}20\rm{\,kpc}$ for
${\sigma_{DM}}/{m_{X}}=1\rm{\,cm^{2}}{gr^{-1}}$, in accord with the range of
the spatial separations shown in the left panels of Figure 11 for the
${\sigma_{DM}}/{m_{X}}=1\rm{\,cm^{2}}{gr^{-1}}$ simulations. As noticed by the
authors, these bounds are an order of magnitude smaller than the measured
offsets for these merging systems ($\sim 100-300\rm{\,kpc}$), and thus in
tension with their simulations. Their conclusions were however affected by the
large observational uncertainties for the measured offsets in massive mergers
of $\sim 10^{15}\rm{\,M_{\odot}}$.

Nonetheless, according to Kim et al. (2021) the presence of offsets between
the different mass components of El Gordo cluster is statistical significant,
which in turn implies that the null hypothesis of zero size offsets is
statistically disfavored.

On the other hand is however difficult to assess the statistical significance
of the offset sizes we found for the merger model XDBf_sb of Figure 12 with
${\sigma_{DM}}/{m_{X}}=4\rm{\,cm^{2}}{gr^{-1}}$, this indeterminacy being
caused by the absence of positional errors for the various peak locations used
to calculate the offsets.

In order to minimize the impact of statistical uncertainties we now attempt to
make a rough estimate of the positional errors for the offsets having the
largest sizes. From the relative location of the various peaks, as extracted
from Figure 6 of Kim et al. (2021), it is easily seen that the galaxy number
density peak of the NW cluster is spatially offset from its parent mass
centroid by about $\sim 120\rm{\,kpc}$. However, for a galaxy group the
determination of the galaxy number density peak location is heavily affected
by shot noise. For a group with $\sim 100$ member galaxies, Kim et al. (2017)
extrapolate to $\sim 160\rm{\,kpc}$ the expected uncertainty.

The large offset ($d^{NW}_{SZ-DM}\sim 150\rm{\,kpc}$) between the SZ emission
peak and the NW mass centroid is similarly influenced by the large error in
the SZ centroid estimate ($\sigma_{SZ}\sim 70\rm{\,kpc}\sim
1^{\prime}.4/(S/N)$, ZYL15). Moreover, for the SE cluster the BCG-DM offset
($\sim 60\rm{\,kpc}$) is approximately within the $\sim 1\sigma$ range of the
mass centroid uncertainty: $\sigma_{DM}\sim 40\rm{\,kpc}$. We have crudely
estimated $\sigma_{DM}$ from Kim et al. (2021), according to the authors a
$2\sigma$ uncertainty in the position of the mass centroids approximately
corresponds to $\sim 10^{{}^{\prime\prime}}$ ($\sim 80\rm{\,kpc}$ at
$z=0.87$).

The size of positional errors for the X-ray emission peak of the SE cluster is
expected to be relatively less significant than for the other centroids, due
to the squared dependency of the X-ray emission with the gas density. These
errors are quantified according to the observational procedure adopted to
process the raw X-ray images and to extract the positions of substructures in
the X-ray emission, together with the relative errors. In the case of the El
Gordo cluster studied here, the effective positional error of the X-ray peak
is relatively small and set by the angular resolution ($\sim
0.5^{{}^{\prime\prime}}$) of the Chandra X-ray image 222J. Kim: private
communication.

The uncertainty in the size of the observed offset between the X-ray and the
SE mass peak is then dominated the uncertainty $\sigma_{DM}$ in the position
of the mass peak, this allows us to put the offset in the range
$d^{SE}_{X-DM}\sim 100\pm 40\rm{\,kpc}$. An important result which follows
from this finding is that SIDM merger simulations of the El Gordo cluster with
${\sigma_{DM}}/{m_{X}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2\rm{\,cm^{2}}{gr^{-1}}$ are
marginally inconsistent with the observed offset $d^{SE}_{X-DM}$. As it can be
seen from the maps of Figure 11, for these models the offset of the X-ray peak
is always limited to $d^{SE}_{X-DM}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}40\rm{\,kpc}$.

In the following we will thus adopt the view that the estimated offset sizes
for the El Gordo cluster are in best accord with SIDM merger simulations for
values of the DM cross-section in the range ${\sigma_{DM}}/{m_{X}}\sim
4-5\rm{\,cm^{2}}{gr^{-1}}$.

To resolve the inconsistency between this range of values for
${\sigma_{DM}}/{m_{X}}$ and its current upper limits on cluster scales
(${\sigma_{DM}}/{m_{X}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1\rm{\,cm^{2}}{gr^{-1}}$), we
suggest here that the SIDM model proposed so far should be considered as only
a first order approximation to the modeling of the physical processes
describing DM interaction. Specifically, we argue that the capacity of DM to
exhibit collisional properties during a cluster merging is correlated with
some other feature of the collision itself.

A simple logical choice to explain this behavior of DM during the mergings is
that the DM collisional properties are not always present, but are switched on
in accordance with the energy of the cluster collisions. This clearly implies
that an SIDM model based on scattering between DM particles should be regarded
as a sort of low order approximation to a much more complex physical process,
in which DM interactions between the two colliding DM halos come into play
according to some energy threshold $E_{\rm crit}$. For a merging cluster with
a collisional energy below this threshold, then the DM of the two halos should
stay in a sort of ground state and during the cluster collision will exhibit
the usual collisionless properties.

It must be stressed such an energy behavior is completely different from the
mechanism previously proposed (Kaplinghat et al. 2016) , in which a power-law
dependency of $\sigma_{DM}/m_{X}$ on the collision velocity has been
postulated in order to satisfy its present upper limits over a wide range of
scales.

To provide a viable theoretical framework for the underlying physical
processes leading to such a step-like behavior for DM is beyond the scope of
this work, nonetheless we can try to assess the validity of this hypothesis by
looking at whether there is some dependency of the observed offsets on the
energy of the merging system.

To validate the proposed picture we take as reference the El Gordo cluster
studied here, for which the measured galaxy-DM and BCG-DM offsets lie in the
range $50\sim 100\rm{\,kpc}$. The best match to the data is obtained for
models XDBf_sa and XDBf_sb, that have initial collision parameters
$\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
V,\leavevmode\nobreak\ d_{ini}\\}$ of model Bf_rc20 : $\\{1.6\cdot
10^{15}\rm{\,M_{\odot}},\leavevmode\nobreak\ 2.3,\leavevmode\nobreak\
2,500\rm{\,km\leavevmode\nobreak\ s^{-1}},\leavevmode\nobreak\
6.1\rm{\,Mpc}\\}$, where $d_{ini}=2(r^{1}_{200}+r^{2}_{200})\simeq
2(1.74+1.32)$. With these parameters we can estimate from Equation (17) the
energy of the collision: $E_{\rm{EG}}\sim 1.4\cdot 10^{64}{\rm{\,ergs}}$.

We now contrast this finding against a merging system for which no appreciable
offsets have been detected between the DM and the stellar component,
specifically we will use the Bullet Cluster. This is a well-known cluster
merger which has been the subject of many numerical studies (Springel & Farrar
2007; Mastropietro & Burkert 2008; Lage & Farrar 2014, 2015; Robertson et al.
2017). We extract the collision parameters
$\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
V,\leavevmode\nobreak\ d_{ini}\\}$ from Table 1 (LF entry) of Lage & Farrar
(2015): $\\{1.9\cdot 10^{15}\rm{\,M_{\odot}},\leavevmode\nobreak\
7.3,\leavevmode\nobreak\ 2,800\rm{\,km\leavevmode\nobreak\
s^{-1}},\leavevmode\nobreak\ 2.8\rm{\,Mpc}\\}$. According to these collision
parameters, the total energy content (17) of the Bullet Cluster is
$E_{\rm{Bullet}}\sim 3\cdot 10^{63}{\rm{\,ergs}}$.

According to our hypothesis, these simple estimates bracket the allowed values
of $E_{\rm crit}$ to be approximately in the range $E_{\rm{Bullet}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}E_{\rm crit}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}E_{\rm{EG}}$. To further
corroborate these findings, it would be interesting to verify if there are
other major merger clusters that exhibit large galaxy-DM peak offsets ($\sim
100\rm{\,kpc}$). If the merging clusters are as massive as El Gordo, this
would clearly lend support to the idea that there is a critical energy
threshold $E_{\rm crit}$ around the value of $E_{\rm{EG}}$.

From Table 2 of Kim et al. (2017), it appears that the Sausage Cluster CIZA
J2242.8+5301 at $z=0.19$ is the only merging cluster which fits these
constraints, apart from El Gordo. This cluster is characterized by an
elongated X-ray emission along the north-south merger axis, hence the name.
The morphology is suggestive of an almost head-on collision, with
spectroscopic data (Dawson et al. 2015) favoring a collision plane close to
that of the sky. The total estimated mass (Jee et al. 2015) is of about $\sim
2\cdot 10^{15}\rm{\,M_{\odot}}$, with a mass ratio close to unity and the DM
centroids separated by about $\sim 1\rm{\,Mpc}$.

For this cluster Table 2 of Kim et al. (2017) reports galaxy-DM offsets
between $\sim 50\rm{\,kpc}$ and $\sim 300\rm{\,kpc}$, although with large
uncertainties. In the context of SIDM it should however be considered
particularly indicative that for the northern group the galaxy centroid
appears leading the corresponding mass peak. According to Jee et al. (2015),
the detection of this offset is highly significant, while for the southern
cluster the errors are relatively large.

We estimate the collision parameters of the Sausage Cluster from Jee et al.
(2015): $\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
V,\leavevmode\nobreak\ d_{ini}\\}=\\{1.1\cdot
10^{15}\rm{\,M_{\odot}},\leavevmode\nobreak\ 1.12,\leavevmode\nobreak\
1,300\rm{\,km\leavevmode\nobreak\ s^{-1}},\leavevmode\nobreak\
4\rm{\,Mpc}\\}$, where $d_{ini}=2(r^{1}_{200}+r^{2}_{200})\sim 4\rm{\,Mpc}$
and $V$ is estimated from Equation (14) of Jee et al. (2015) by setting $b=0$
and $t_{\rm impact}\sim 11$ Gyr. We can now find the energy of the merging,
which turns out to be about $E_{\rm{Sausage}}\sim 1.5\cdot
10^{64}{\rm{\,ergs}}$.

This energy value is very close to that of $E_{\rm{EG}}$ and strongly supports
the hypothesis of the existence of an energy threshold which regulates DM
behavior, thus implying that during the cluster collision DM can exhibit
collisional properties according to the amount of collisional energy.

Numerical simulations of the Sausage Cluster have been previously performed by
some authors (van Weeren et al. 2011; Donnert et al. 2017; Molnar & Broadhurst
2017), but without including the effect of SIDM. From the SIDM merging
simulations of Kim et al. (2017) the expected galaxy-DM offset for merging
clusters of $\sim 10^{15}\rm{\,M_{\odot}}$ is of the order of $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}20\rm{\,kpc}$ for
${\sigma_{DM}}/{m_{X}}=1\rm{\,cm^{2}}{gr^{-1}}$, which is smaller than the
estimated offset for the Sausage cluster by about a factor of ten.

We argue that these inconsistencies would be solved in the SIDM scenario
proposed here, and that the study of the different offsets in major cluster
mergers is going to provide significant clues for the nature of DM. In this
context SIDM merging simulations of the Sausage Cluster, in connection with
the results extracted here from the SIDM runs of the El Gordo cluster, will
most likely lead to meaningful insights into the modeling of DM self-
interactions based on particle collisions.

###### Acknowledgements.

The author is grateful to Congyao Zhang for his many clarifying comments about
his paper that have helped to produce this work. The author would also like to
thank Jinhyub Kim for helping with Figure 12, John Miller is also acknowledged
for carefully reading and improving the manuscript. The computations of this
paper were carried out using the Ulisse cluster at SISSA and the Marconi
cluster at CINECA (Italy), under a SISSA-CINECA agreement.

##  Appendix A Parallel implementation of the DM scattering algorithm and code
validation

We now describe the parallelization of the Monte Carlo scheme described in
Section 2.3. Let us consider the whole set $\\{i\\}{\subseteq\Omega}$ of the
$N_{DM}$ DM particles which are distributed across the whole computational
volume $\Omega$. In a serial code scattering events will be sampled by
applying Equation (19) to the whole set $\\{i\\}{\subseteq\Omega}$ of DM
particles. However the scattering process is cumulative, so that any
scattering event requires knowledge from previous events. If one implements
the algorithm in a serial code this condition is implicit, but with a parallel
code this is not obvious.

The Monte Carlo method is then implemented as follows. Let us assume that the
parallel code subdivides the computation among $P$ tasks, and the simulation
volume is correspondingly split into $\Omega_{p}$ non-overlapping subvolumes,
with $p=0,1,\ldots,P-1$ being the task index. The load balancing domain
decomposition is done by cutting segments of a globally ordered Hilbert space-
filling curve, as in the code GADGET-2 (Springel 2005).

We now evaluate the scattering probability (19) for the subset
$\\{i\\}{\subseteq\Omega_{p}}$ of DM particles which are members of the task
$p$. If a scattering event occurs then the tagged neighbor particle $j$ could
be either a member of $p$ or of another task $q$, in the latter case the
particle having been imported by $p$ during the root-finding neighbor search
(18). Moreover, even if $j$ is a member of the local task $p$, it could still
be a neighbor member of a particle in another task $q$. The same applies to
particle $i$. Therefore, a blind implementation of the Monte Carlo algorithm
in a parallel code could lead to the same particle being scattered two or more
times with the same initial velocity.

To avoid this problem we then proceed as follows. Let us first assume that at
any given time all the scattering pairs $\\{i,j\\}$ of the $P$ tasks have been
identified. For the generic pair $\\{i,j\\}_{p}$ in the task $p$ there is then
the possibility that the neighbor particle $j$ is a member of another task
$q$, otherwise either $i$ or $j$ or both could be neighbors of another
scattering pair in the task $q$. We now denote as $\\{i,j\\}^{in}_{p}$ all of
the scattering pairs of $p$ which do not satisfy any of the above conditions,
i.e. the scattering takes place entirely within the task $p$. We then split
the set of scattering pairs $\\{i,j\\}_{p}$ in the task $p$ between the subset
$\\{i,j\\}^{in}_{p}$ and all the other pairs, which we denote as
$\\{i,j\\}^{oth}_{p}$.

To avoid multiple scatterings it is now necessary that each task $p$ knows if
a member of a scattering pair $\\{i,j\\}^{oth}_{p}$ has undergone a scattering
event in another task, this in order to accordingly update its velocity. To
accomplish this for each pair $\\{i,j\\}^{oth}_{p}$ we associate a 64 bit
integer key defined as $ID^{(i_{K},j_{K})}_{p}=i_{K}(i)+Mj_{K}(j)$, where
$i_{K}(i)$ denotes an integer key of the local DM particle $i$. The integer
key $i_{K}(i)$ is unique for each DM particle and it is assigned initially, it
is a DM particle property which moves with the particle when this migrates
between processors. The integer $M$ is is a large integer ($2^{31}$) such that
$M>>N_{DM}$, so that in the $P$ tasks the integers $ID^{(i_{K},j_{K})}_{p}$
are unique for each scattering pair.

To consistently incorporate multiple scatterings in the parallel code, we now
perform a global sorting of the $P$ sets of integers $ID^{(i_{K},j_{K})}_{p}$,
with $p=0,1,\ldots,P-1$. After the sorting each task $p$ knows the global list
of sorted keys, and therefore all of the other scatterings and subsequently in
which task they are taking place.

Finally, the post-scattering velocities (20) of each scattering pair in the
task $p$ are evaluated as follows. After the global sorting each task will
process its local set of scattering pairs by looping iteratively over its
local list of sorted keys, this in increasing key order because of the
sorting. Pairs which exhibit a scattering dependency with
$ID^{(i_{K},j_{K})}_{p}>ID^{(j_{K},m_{K})}_{q}$ are put in a waiting list
until the updated velocities of particle $j$ have been received from task $q$.
Note that this implementation allows for a single particle to undergo multiple
scatterings.

To improve code performances, we make use of asynchronous communications to
optimize data transfer between processors. We exploit the advantage of using
this communication mode by splitting the list of local pairs
$\\{i,j\\}^{in}_{p}$ and by scattering a subset of these pairs whenever the
task is waiting updated particle velocities from other tasks.

To test the SIDM part of the code we evolve an isolated DM NFW halo, initially
at equilibrium, for different values of the self-interacting DM cross section.
We have considered the following values:
${\sigma_{DM}}/{m_{X}}=0,\leavevmode\nobreak\ 1,$ and
$10\rm{\,cm^{2}}{gr^{-1}}$. For the DM halo we set $M_{200}=1.3\cdot
10^{15}\rm{\,M_{\odot}}$ at $z=0.87$, so that $R_{200}=1.63$ Mpc and for the
concentration parameter $c_{200}=2.55$ from Equation (5). The initial
condition setup of the DM halo is described in Section 2.2.1.

For the different runs we show here at a given time $t$ the spherically
averaged DM radial density profile $\rho_{DM}(r,t)$ and the scattering rate
per particle $\Gamma(r,t)$. The theoretical expression of $\Gamma(r,t)$ at the
time $t$ is evaluated as (Robertson et al. 2015)

$\Gamma_{th}(r)=\frac{\sigma_{DM}}{m_{X}}<v_{pair}(r)>\rho_{DM}(r)$ (28)

where $\rho_{DM}(r)$ is the local DM density and $<v_{pair}(r)>$ the mean
pairwise velocity. Assuming a Maxwell-Boltzmann distribution this quantity can
be written as (Robertson et al. 2015)
$<v_{pair}(r)>=(4/\sqrt{\pi})\sigma_{1D}(r)$, where $\sigma_{1D}(r)$ is the
one-dimensional velocity dispersion of the NFW halo (Lokas & Mamon 2001).

To contrast the predicted rate $\Gamma_{th}(r)$ with the corresponding rate
$\Gamma_{sim}(r)$ extracted from the simulations, we measure $\Gamma_{sim}(r)$
as follows. We first accumulate in radial bins the number of scattering events
as a function of radius, together with the number of particles in the same
bin. At the epoch $t$ the binned scattering rate $\Gamma_{sim}(r)$ is then
obtained by dividing these two quantities. It must be stressed that during the
simulations the rates $\Gamma_{sim}(r)$ are obtained by evaluating the local
scattering probability (19), but without updating the particle velocities (20)
of the pair. This is to avoid any modification to the DM density profile,
which otherwise will make Equation (28) not applicable.

The results we obtain after $t=1$ Gyr are shown in Figure 13, in which for the
different cases are plotted the DM density profiles (left panel) and the
corresponding scattering rates (right panel). For completeness, in the left
panel we also show the density profile (solid black line) of a halo without
SIDM (${\sigma_{DM}}/{m_{X}}=0$). As expected, there is the development of
cored density profiles for which at inner radii the flattening of density is
higher as ${\sigma_{DM}}/{m_{X}}$ increases. For the considered SIDM cases,
the scattering rates are in full agreement with the theoretical predictions.

The results of similar tests have been shown by various authors (Vogelsberger
et al. 2012; Robertson et al. 2015, 2017; ZuHone et al. 2019; Fisher et al.
2021), we compare our density profile results with those of a similar test
carried out by ZuHone et al. (2019). The authors adopt the same SIDM
implementation (Vogelsberger et al. 2012) employed here, and moreover assume
the same approximations for the DM cross-section ${\sigma_{DM}}$. In their
Section 3.1 the authors investigate the time evolution of an isolated halo
comprising gas and DM. They assume a cluster gas fraction of $0.12$, we then
consider dynamically negligible the impact of the gas when comparing at inner
radii the different DM density profiles.

For the halo DM density profile ZuHone et al. (2019) chose an Hernquist
profile with $M^{Z}_{200}=1.25\cdot 10^{15}\rm{\,M_{\odot}}$ and
$a^{Z}_{H}=600\rm{\,kpc}$ for the cluster mass and profile scale radius,
respectively. Our halo mass approximately coincides with $M^{Z}_{200}$, so
that our rescaled profiles can be directly compared with the corresponding
ones of ZuHone et al. (2019).

Following Springel (2005) we relate the scale radius $r_{s}$ and concentration
parameter $c$ of the NFW profile to the scale length $a_{H}$ of an Hernquist
density profile with the same mass within $r_{200}$

$a_{H}=r_{s}\sqrt{2[\ln(1+c)-c/(1+c)]}.$ (29)

For the halo parameters we use $r_{s}=r_{200}/c\simeq 0.64\rm{\,Mpc}$ and
$a_{H}\simeq 660\rm{\,kpc}$, which is only $\sim 10\%$ higher than
$a^{Z}_{H}$. Thus at inner radii our profiles can be compared with the DM halo
profiles displayed in top left panel of Figure 1 of ZuHone et al. (2019). From
the left panel of Figure 13 at the initial epoch ($t=0$) we estimate at the
radius $r\sim 10^{-2}r_{200}\simeq 16\rm{\,kpc}$ a central DM density of
$\rho_{DM}/\rho_{c}\simeq 8\cdot 10^{4}$. For the chosen redshift ($z=0.87$)
this value corresponds to $\rho_{DM}\simeq 3\cdot
10^{7}\rm{\,M_{\odot}}kpc^{-3}$, in accord with the corresponding value
estimated from Figure 1 of ZuHone et al. (2019).

We now estimate at $t=1$ Gyr the central value of the halo DM density for the
SIDM run with ${\sigma_{DM}}/{m_{X}}=10\rm{\,cm^{2}}{gr^{-1}}$. For this test
from Figure 13 we obtain at the radius $r\sim 10^{-2}r_{200}$ a central DM
density of $\rho_{DM}/\rho_{c}\simeq 1.8\cdot 10^{4}$, which is equivalent to
$\rho_{DM}\simeq 6.6\cdot 10^{6}\rm{\,M_{\odot}}kpc^{-3}$. For the same SIDM
run this value is in close agreement with the corresponding DM density,
estimated at the same radius and at the same epoch from Figure 1 of ZuHone et
al. (2019). The results of our tests then validate our implementation of DM
self-interactions, and demonstrate that the resulting code can be safely used
to perform the SIDM simulations presented here.

Figure 13:  Radial density (left panel) and scattering rate (right panel)
profiles are shown at $t=1$ Gyr for three SIDM runs with different DM
scattering cross-section: ${\sigma_{DM}}/{m_{X}}=0,\leavevmode\nobreak\ 1,$
and $10\rm{\,cm^{2}}{gr^{-1}}$. The tests were carried out for a NFW DM halo
with $M_{200}=1.3\cdot 10^{15}\rm{\,M_{\odot}}$ and $R_{200}=1.63$ Mpc,
initially at equilibrium. In the left panel we also show the density profile
of a run without SIDM, ${\sigma_{DM}}/{m_{X}}=0$, at $t=0$ ( dash black line)
and at $t=1$ Gyr ( solid black line). The solid blue line is the NFW
theoretical profile.

## References

  * Allen et al. (2008) Allen, M. G., Groves, B. A., Dopita, M. A. et al., 2008, ApJS, 178, 20
  * Agertz et al. (2007) Agertz, O., Moore, B., Stadel, J. et al., 2007, MNRAS, 380, 963
  * Anders & Grevesse (1989) Anders, E. & Grevesse, N., 1989 , Geochim. Cosmochim. Acta, 53, 197
  * Ascasibar & Markevitch (2006) Ascasibar, Y. & Markevitch, M., 2006, ApJ, 650, 102
  * Asencio et al. (2021) Asencio, E., Banik, I. & Kroupa, P., 2021, MNRAS, 500, 5249
  * Asencio et al. (2023) Asencio, E., Banik, I. & Kroupa, P., 2023, arXiv:2308.00744
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S., 1987, ”Galactic dynamics”, Princeton Univeristy Press
  * Buote (2002) Buote, D. A., 2002, Astrophysics and Space Science Library, 272 , 79, Feretti, L., Gioia, I. M. and Giovannini, G. (eds), Kluwer Academic Publishers
  * Burkert (1995) Burkert, A., 1995, ApJ, 447, L25
  * Canning et al. (2012) Canning, R. E. A., Russell, H. R., Hatch, N. A. et al., 2012, MNRAS, 420, 2956
  * Cavaliere & Fusco-Femiano (1978) Cavaliere, A. & Fusco-Femiano, R., 1978, A&A, 70, 677
  * Chadayammuri et al. (2022) Chadayammuri, U., ZuHone, J., Nulsen, P. et al., 2022, MNRAS, 509, 1201
  * Clowe et al. (2006) Clowe, D., Bradač, M., Gonzalez, A. H. et al., 2006, ApJ, 648, L109
  * Cross et al. (2023) Cross, D., Thoron, G., Jeltema, T. E., et al., 2023, arXiv:2304.10128
  * Dawson et al. (2015) Dawson, W. A., Jee, M. J., Stroe, A. et al., 2015, ApJ, 805 143
  * Dehnen & Aly (2012) Dehnen, W. & Aly, H., 2012, MNRAS, 425, 1068
  * Depropis et al. (2021) De Propris, R., West, M. J., Andrade-Santos, F., et al., 2021, MNRAS, 500, 310
  * Diego et al. (2020) Diego, J. M., Molnar, S. M., Cerny, C. et al., 2020, ApJ, 904, 106
  * Donnert (2014) Donnert, J. M. F., 2014, MNRAS, 438, 1971
  * Donnert et al. (2017) Donnert, J. M. F., Beck, A. M., Dolag, K. & Röttgering, H. J. A., 2017, MNRAS, 471, 4587
  * Duffy et al. (2008) Duffy, A. R., Schaye, J., Kay, S. T. & Dalla Vecchia, C., 2008, MNRAS, 390, L64-L68
  * Fisher et al. (2021) Fischer, M. S., Brüggen, M., Schmidt-Hoberg, K. et al., 2021, MNRAS, 505, 851
  * Fisher et al. (2022) Fischer, M. S., Brüggen, M., Schmidt-Hoberg, K. et al., 2022, MNRAS, 510, 4080
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q., MacNeice, P., Rosner, R., Truran, J. W. & Tufo, H., 2000 ApJS, 131, 273
  * García-Senz et al. (2012) García-Senz, D., Cabezón, R. M. & Escartín, J. A., 2012, A&A, 538, A9
  * Golovich et al. (2017) Golovich, N., van Weeren, R. J., Dawson, W. A., Jee, M. J. & Wittman, D., 2017, ApJ, 838, 110
  * Halbesma et al. (2019) Halbesma, T. L. R., Donnert, J. M. F., de Vries, M. N. & Wise, M. W., 2019, MNRAS, 483, 3851
  * Hallman et al. (2004) Hallman, E. J. & Markevitch, M. , 2004, ApJ, 610, L81
  * Harvey et al. (2014) Harvey, D., Titley, E., Massey, R., et al., 2014, MNRAS, 441, 404
  * Harvey et al. (2019) Harvey, D., Robertson, A., Massey, R. & McCarthy, I. G., 2019, MNRAS, 488, 1572
  * Hasselfield et al. (2013) Hasselfield, M., Hilton, M. & Marriage, T., A., 2013, J. Cosmology Astropart. Phys., 7, 8
  * Heinz et al. (2003) Heinz, S., Churazov, E., Forman, W., Jones, C. & Briel, U. G., 2003, MNRAS, 346, 13
  * Hernquist & Katz (1989) Hernquist, L. & Katz, N., 1989 ApJS, 70, 419
  * Hockney & Eastwood (1988) Hockney, R. W. & Eastwood, J. W., 1988, ”Computer simulation using particles”, Bristol: Hilger
  * Itoh et al. (1998) Itoh, N., Kohyama, Y. & Nozawa, S., 1998, ApJ, 502, 7
  * Jee et al. (2014) Jee, M. J., Hughes, J. P., Menanteau, F., et al., 2014, ApJ, 785, 20
  * Jee et al. (2015) Jee, M. J. Stroe, A., Dawson, W., et al., 2015, ApJ, 802, 46
  * Kahlhoefer et al. (2014) Kahlhoefer, F., Schmidt-Hoberg, K., Frandsen, M. T. & Sarkar, S., 2014, MNRAS, 437, 2865
  * Kaplinghat et al. (2016) Kaplinghat, M., Tulin, S. & Yu, H.-B., 2016, Phys. Rev. Lett., 116, 041302
  * Kazantzidis et al. (2004) Kazantzidis, S., Magorrian, J. & Moore, B., 2004, ApJ, 601, 37
  * Kim et al. (2017) Kim, S. Y., Peter, A. H. G. & Wittman, D. et al., 2017, MNRAS, 469, 1414
  * Kim et al. (2021) Kim, J., Jee, M. J., Hughes, J. P. et al. , 2021, ApJ, 923, 101
  * Kravtsov et al. (2018) Kravtsov, A. V., Vikhlinin, A. A. &Meshcheryakov, A. V., 2018, Astronomy Letters, 44, 8
  * Lage & Farrar (2014) Lage, C. & Farrar, G. R., 2014, ApJ, 787, 144
  * Lage & Farrar (2015) Lage, C. & Farrar, G. R., 2015, J. Cosmology Astropart. Phys., 2015, 38
  * Lee & Komatsu (2010) Lee, J. & Komatsu, E., 2010, ApJ, 718, 60
  * Lindnet et al. (2014) Lindner, R. R., Baker, A. J., Hughes, J. P. et al, 2014, ApJ, 786, 49
  * Lokas & Mamon (2001) Łokas, E. L. & Mamon, G. A., 2001, MNRAS, 321, 155
  * Machado & Lima Neto (2013) Machado, R. E. G. & Lima Neto, G. B., 2013, MNRAS, 430, 3249
  * Machado et al. (2015) Machado, R. E. G., Monteiro-Oliveira, R., Lima Neto, G. B. & Cypriano, E. S., 2015, MNRAS, 451, 3309
  * Markevitch et al. (2004) author = Markevitch, M., Gonzalez, A. H., Clowe, D. et al., 2004, ApJ, 606, 819
  * Martel et al. (2014) Martel, H., Robichaud, F. & Barai, P., 2014, ApJ, 786, 79
  * Mastropietro & Burkert (2008) Mastropietro, C. & Burkert, A., 2008, MNRAS, 389, 967
  * Mazzotta et al. (2004) Mazzotta, P., Rasia, E., Moscardini, L. & Tormen, G., 2004, MNRAS, 354, 10
  * McCarthy et al. (2007) McCarthy, I. G., Bower, R. G., Balogh, M. L., et al., 2007, MNRAS, 376, 497
  * Menanteau et al. (2010) Menanteau, F., González, J., Juin, J.-B. et al., 2010, ApJ, 723, 1523
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J. & Terzić, B., 2006, AJ, 132, 2685
  * Molnar et al. (2013) Molnar, S. M., Chiu, I.-N. T., Broadhurst, T. & Stadel, J. G., 2013, ApJ, 779, 63
  * Molnar & Broadhurst (2015) Molnar, S. M. &Broadhurst, T., 2015, ApJ, 800, 37
  * Molnar & Broadhurst (2017) Molnar, S. M. &Broadhurst, T., 2017, ApJ, 841, 46
  * Molnar & Broadhurst (2018) Molnar, S. M. &Broadhurst, T., 2018, ApJ, 862, 112
  * Moura et al. (2021) Moura, M. T., Machado, R. E. G. & Monteiro-Oliveira, R., 2021, MNRAS, 500, 1858
  * Markevitch et al. (2002) Markevitch, M., Gonzalez, A. H., David, L., et al., 2002, ApJ, 567, L27
  * Markevitch & Vikhlinin (2007) Markevitch , M. & Vikhlinin, A., 2007 , Phys. Rep., 443, 1
  * Menanteau et al. (2012) Menanteau, F., Hughes, J. P., Sifón, C. et al., 2012, ApJ, 748, 7
  * Mitchell et al. (2009) Mitchell, N. L., McCarthy, I. G., Bower, R. G., Theuns, T. & Crain, R. A., 2009 , MNRAS, 395, 180
  * Molnar et al. (2012) Molnar, S. M., Hearn, N. C. & Stadel, J. G., 2012, ApJ, 748, 45
  * Molnar (2016) Molnar, S., 2016 , Front. Astr. Space Sci., 2, 7
  * Nagai et al. (2007) Nagai, D., Kravtsov, A. V. & Vikhlinin, A., 2007, ApJ, 668, 1
  * Ng et al. (2015) Ng, K. Y., Dawson, W. A., Wittman, D. et al., 2015, MNRAS, 453, 1531
  * Poole et al. (2006) Poole, G. B., Fardal, M. A., Babul, A., et al, 2006, MNRAS, 373, 881
  * Poole et al. (2008) Poole, G. B., Babul, A., McCarthy, I. G., Sanderson, A. J. R. & Fardal, M. A., 2008, MNRAS, 391, 1163
  * Power et al. (2003) Power, C., Navarro, J. F., Jenkins, A. et al, 2003, MNRAS, 338, 14
  * Price (2012) Price, D. J., 2012, J. Comp. Phys., 231, 759
  * Ragozzine et al. (2012) Ragozzine, B., Clowe, D., Markevitch, M., Gonzalez, A. H. & Bradač, M., 2012, ApJ, 744, 94
  * Randall et al. (2008) Randall, S. W., Markevitch, M, Clowe, D., Gonzalez, A. H. & Bradač, M., 2008, ApJ, 679, 1173
  * Ricker (1998) Ricker, P. M., 1998, ApJ, 496, 670
  * Ricker & Sarazin (2001) Ricker, P. M. & Sarazin, C. L., 2001, ApJ, 561, 621
  * Ritchie & Thomas (2002) Ritchie, B. W. & Thomas, P. A., 2002, MNRAS, 329, 675
  * Robertson et al. (2015) Robertson, A., Massey, R., Eke, V. & Bower, R., 2015, MNRAS, 453, 2267
  * Robertson et al. (2017) Robertson, A., Massey, R. & Eke, V., 2017, MNRAS, 465, 569
  * Robertson et al. (2019) Robertson, A., Harvey, D., Massey, R. et al., 2019, MNRAS, 488, 3646
  * Roediger et al. (2012) Roediger, E., Lovisari, L., Dupke, R. et al, 2012, MNRAS, 420, 3632
  * Roediger et al. (2013) Roediger, E., Kraft, R. P., Forman, W. R., Nulsen, P. E. J. & Churazov, E., 2013, ApJ, 764, 60
  * Roettiger et al. (1996) Roettiger, K., Burns, J. O. & Loken, C., 1996, ApJ, 473, 651
  * Russell et al. (2010) Russell, H. R., Sanders, J. S., Fabian, A. C., et al., 2010, MNRAS, 406, 1721
  * Ruszkowski & Oh (2010) Ruszkowski, M. & Oh, S. P., 2010, ApJ, 713, 1332
  * Sarazin (2002) Sarazin, C. L., 2002, Astrophysics and Space Science Library, 272 , 1, Feretti, L., Gioia, I. M. and Giovannini, G. (eds), Kluwer Academic Publishers
  * Schmidt et al. (2017) Schmidt, W., Byrohl, C., Engels, J. F., Behrens, C. & Niemeyer, J. C., 2017, MNRAS, 470, 142
  * Sheardown et al. (2018) Sheardown, A., Roediger, E., Su, Y. et al., 2018, ApJ, 865, 118
  * Sheardown et al. (2019) Sheardown, A., Fish, T. M., Roediger, E., et al., 2019, ApJ, 874, 112
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M. et al., 2003, MNRAS, 343, 978
  * Shen et al. (2022) Shen, X., Brinckmann, T., Rapetti, D. et al., 2022, MNRAS, 516, 1302
  * Springel (2005) Springel, V., 2005, MNRAS, 364, 1105
  * Springel (2005) Springel, V., Di Matteo, T. & Hernquist, L. 2005, MNRAS, 361, 776
  * Springel & Farrar (2007) Springel, V. & Farrar, G. R., 2007, MNRAS, 380, 911
  * Tulin & Yu (2018) Tulin, S. & Yu, H.-B., 2018, Phys. Rep., 730, 1
  * Valdarnini (2006) Valdarnini, R., 2006, New Astronomy, 12, 71
  * Valdarnini (2016) Valdarnini, R., 2016, ApJ, 831, 103
  * Valdarnini (2019) Valdarnini, R., 2019, ApJ, 874, 42
  * Valdarnini & Sarazin (2021) Valdarnini, R. & Sarazin, C. L., 2021, MNRAS, 504, 5409
  * van den Bosch et al. (2005) van den Bosch, F. C., Weinmann, S. M., Yang, X., 2005, MNRAS, 361, 1203
  * Vogelsberger et al. (2012) Vogelsberger, M., Zavala, J. & Loeb, A., 2012, MNRAS, 423, 3740
  * Voit (2005) Voit, G.M. 2005, RMP, 77, 207
  * van Weeren et al. (2011) van Weeren, R. J., Brüggen, M., Röttgering, H. J. A. & Hoeft, M., 2011, MNRAS, 418, 230
  * Wittman et al. (2018) Wittman, D., Golovich, N. & Dawson, W. A., ApJ, 869, 10
  * Zhang et al. (2014) Zhang, C., Yu, Q. & Lu, Y., 2014, ApJ, 796, 138
  * Zhang et al. (2015) Zhang, C., Yu, Q. & Lu, Y., 2015, ApJ, 813, 129 (ZYL15)
  * Zhang et al. (2018) Zhang, C., Yu, Q. & Lu, Y., 2018, ApJ, 855, 36
  * Zhao et al. (2004) Zhao, P., Jerius, D. H., Edgar, R. J. et al., 2004, Proc SPIE, 5165, 482
  * Zitrin et al. (2012) Zitrin, A., Bartelmann, M., Umetsu, K., Oguri, M. & Broadhurst, T., 2012, MNRAS, 426, 2944
  * Zitrin et al. (2013) Zitrin, A., Menanteau, F., Hughes, J. P. et al., 2013, ApJ, 770, L15
  * ZuHone et al. (2009) ZuHone, J. A., Ricker, P. M., Lamb, D. Q. & Karen Yang, H. -Y., 2009, ApJ, 699, 1004
  * ZuHone et al. (2010) ZuHone, J. A., Markevitch, M. & Johnson, R. E., 2010, ApJ, 717, 908
  * ZuHone (2011) ZuHone, J. A., 2011, ApJ, 728, 54
  * ZuHone et al. (2013) ZuHone, J. A., Markevitch, M., Ruszkowski, M. & Lee, D., 2013, ApJ, 762, 69
  * ZuHone et al. (2015) ZuHone, J. A., Kunz, M. W., Markevitch, M., Stone, J. M. & Biffi, V., 2015, ApJ, 798, 90
  * ZuHone & Roediger (2016) ZuHone, J. A. & Roediger, E., 2016, Journal of Plasma Physics, 82, 535820301
  * ZuHone et al. (2019) ZuHone, J. A., Zavala, J. & Vogelsberger, M., 2019, ApJ, 882, 119