11institutetext: SISSA/ISAS, via Bonomea 265, I-34136 Trieste, Italy
11email<EMAIL_ADDRESS>
#  An N-body/hydrodynamical simulation study of the merging cluster El Gordo:
a compelling case for SIDM ?

R. Valdarnini 11

(Received ; accepted )

We use a large set N-body/hydrodynamical simulations to study the physical
properties of the merging cluster El Gordo. We found that the observed X-ray
structures, along with other data, can be fairly matched by simulations with
collision velocities $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 impact parameters $600\rm{\,kpc}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}P\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}800\rm{\,kpc}$. The mass of the
primary is constrained to be between $\sim 10^{15}\rm{\,M_{\odot}}$ and $\sim
1.6\cdot 10^{15}\rm{\,M_{\odot}}$, in accord with recent lensing-based mass
measurements. Moreover, a returning, post-apocenter, scenario is not supported
by our head-on simulations.

We also consider merger models that incorporate dark matter self-interactions.
The simulation results show that the observed spatial offsets between the
different mass components are well reproduced in self-interacting dark matter
models with an elastic cross-section in the range $\sigma_{DM}/{m_{X}}\sim
4-5\rm{\,cm^{2}}{gr^{-1}}$. In addition, the mean relative line-of-sight
radial velocity between the two brightest cluster galaxies is found to be of
the order of several hundreds of $\rm{\,km\leavevmode\nobreak\ s^{-1}}$. We
argue that these findings provide an unambiguous signature of a dark matter
behavior that exhibits collisional properties in a very energetic high-
redshift cluster collision. The range of allowed values we found for
$\sigma_{DM}/{m_{X}}$ is however inconsistent with present upper limits. To
resolve this tension we suggest the possibility that the self-interacting dark
matter model used here should be considered as only a low order approximation,
and that the underlying physical processes that describe the interaction of
dark matter in major cluster mergers are more complex than can be adequately
represented by the commonly assumed approach based on scattering of dark
matter particles.

###### Key Words.:

— galaxies: clusters: intracluster medium — Hydrodynamics — methods: numerical
— X-rays: galaxies: clusters — cosmology: theory – dark matter

##  1 Introduction

Galaxy clusters are the largest virialized structures known in the Universe
and according to the hierarchical scenario their formation proceeds bottom-up,
through the sequential accretion and merging of smaller structures. They
constitute the most massive virialized objects actually present in the
Universe, with measured virial masses in the range $M\sim
10^{14}-10^{15}\rm{\,M_{\odot}}$. Observations show that about $\sim 80\%$ of
their mass content consists of dark matter (DM), with the remaining $\sim
20\%$ in the form of baryons. The bulk of these baryons ($\sim 90\%$) resides
in the form of an intracluster medium (ICM), with only a few per cent in the
form of galaxies (Voit 2005).

During the process of cluster formation, high-velocity gas flows will
dissipate the infall kinetic energy and the gas will be shock-heated to high
temperatures. Because of the masses involved, at equilibrium these
temperatures will be as high as $T\sim 10^{7}-10^{8}$ K and the ICM will be in
the form of a hot diffuse plasma. Galaxy clusters are thus powerful X-ray
emitters and their X-ray emission can then be used to study the physics of the
ICM and in turn to derive cluster properties with which to constrain the
currently accepted Lambda cold dark matter ($\Lambda$CDM) cosmological model
(Voit 2005).

Moreover observed X-ray maps show that a significant fraction of galaxy
clusters exhibit a disturbed X-ray morphology and are in a dynamically
unrelaxed state (Buote 2002), with two or more subclusters still in the
process of accretion and merging. In this context the merging of galaxy
clusters is essential in a number of ways to probe ICM physics and the
underlying cosmology (Sarazin 2002; Molnar 2016).

Mergers between galaxy clusters can be broadly classified into two categories,
according to the mass ratio between the two main mass subcomponents. Major
mergers occur when the mass ratio of the primary to the secondary is close to
unity, while minor mergers take place whenever the mass of the secondary is
below $\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}30-40\%$ of the
primary mass. X-ray surface brightness maps of the latter reveal the existence
of cold fronts and shocks (Markevitch & Vikhlinin 2007), which can be used to
derive ICM transport properties (Markevitch & Vikhlinin 2007; ZuHone &
Roediger 2016).

In this scenario major mergers of massive galaxy clusters can be considered as
the most energetic events since the Big-Bang, with collision energies in the
range $\sim 10^{63}-10^{64}\rm{\,ergs^{-1}}$ (Sarazin 2002). The energy
involved in these collisions renders these objects unique laboratories to
probe the collisional properties of their components, as well as the
likelihood of such events within the $\Lambda$CDM cosmological model (Lee &
Komatsu 2010; Lage & Farrar 2015; Asencio et al. 2021).

An important effect which will occur during a major merger is the expected
offset between the dissipative ICM and the collisionless components: galaxies
and DM. Additionally, self-interacting DM (SIDM) will also exhibit a further
offset from the galaxy component which will depend on the strength of the
cross-section. Observations of spatial offsets between different components
can then be used to constrain these models (see Tulin & Yu 2018, for a review
).

In recent years the existence of major mergers between galaxy clusters has
been established by a growing number of observations: for example, the Bullet
Cluster 1E0657-56 at redshift $z=0.296$ (Markevitch et al. 2002; Clowe et al.
2006), which demonstrated for the first time the existence of DM by displaying
a clear spatial offset between the collisionless DM component and the X-ray
emitting ICM. Other examples of such systems are the cluster ACT-CL J0102-4915
(Menanteau et al. 2012; Jee et al. 2014; Zitrin et al. 2013), the Sausage
Cluster CIZA J2242.8+5301 (Jee et al. 2015), ZwCl008.8+52 (Golovich et al.
2017), A1758N (Ragozzine et al. 2012) and A2146 (Russell et al. 2010).

Because of the non-linearity of the merging process and the variety of the
associated physical phenomena, N-body/hydrodynamical simulations have proved
to be a very powerful tool for studying and interpreting mergers between
galaxy clusters. In particular, idealized binary merger simulations are
designed to follow, in a Newtonian frame, the collision evolution of two
merging halos which are initially separated and at equilibrium. At variance
with cosmological simulations, this approach allows a single merger to be
studied under carefully controlled initial conditions and can be used to model
a specific merging event (Molnar 2016).

Individual binary cluster merger simulations have been performed by many
authors. These simulations have been aimed either at studying in general the
effects of mergers (Roettiger et al. 1996; Ricker & Sarazin 2001; Ritchie &
Thomas 2002; Poole et al. 2006; McCarthy et al. 2007; Poole et al. 2008;
Mitchell et al. 2009; ZuHone et al. 2009), or at investigating some specific
property of the ICM.

In particular, simulations of merging clusters have been used to study the
origin of cold fronts and gas sloshing (Heinz et al. 2003; Ascasibar &
Markevitch 2006; ZuHone et al. 2010; Roediger et al. 2012; ZuHone & Roediger
2016; Sheardown et al. 2018), the formation of cool cores (ZuHone 2011;
Valdarnini & Sarazin 2021), the offset between X-ray and Sunyaev-Zel’dovich
(SZ) map peaks (Molnar et al. 2012; Zhang et al. 2014) and ICM transport
properties (Ruszkowski & Oh 2010; Roediger et al. 2013; ZuHone et al. 2013,
2015; Schmidt et al. 2017).

Moreover this simulation approach allows simulations to be performed aimed at
studying a particular merging system, so as to contrast the simulation results
with a given observation. Some examples of such specifically designed
simulations include ones of the Bullet Cluster (Springel & Farrar 2007;
Mastropietro & Burkert 2008; Lage & Farrar 2014) and the cluster ACT-CL
J0102-4915 (Donnert 2014; Molnar & Broadhurst 2015; Zhang et al. 2015, 2018).
Other examples are ones of the Sausage cluster (Donnert et al. 2017; Molnar &
Broadhurst 2017), Cygnus A (Halbesma et al. 2019), A2146 (Canning et al. 2012;
Chadayammuri et al. 2022), A3376 (Machado & Lima Neto 2013), A1758N (Machado
et al. 2015), the cluster ZwCl008.8+52 (Molnar & Broadhurst 2018) and A2034
(Moura et al. 2021).

Among these examples of extreme collisions, the cluster ACT-CL J0102-4915 at
$z=0.870$, nicknamed El Gordo, presents at various levels the most interesting
features. It was originally detected as the cluster with the strongest SZ
effect in the Atacama Cosmology Telescope (ACT) survey (Menanteau et al. 2010;
Hasselfield et al. 2013). From optical and X-ray studies Menanteau et al.
(2012) obtained a total mass of the order of $\sim 2\cdot
10^{15}\rm{\,M_{\odot}}$, with a galaxy velocity dispersion $\sigma_{gal}\sim
1,300\rm{\,km\leavevmode\nobreak\ s^{-1}}$ and an integrated temperature in
the range $T_{\mathrm{X}}\simeq 15$ keV. Similar mass estimates have also been
obtained from weak lensing (Jee et al. 2014) and strong lensing (SL) studies
(Zitrin et al. 2013). This range of estimated masses indicates that El Gordo
is a major merger cluster as well as being the most massive cluster at
$z>0.6$, however it was soon realized that the existence of such a massive
cluster at this redshift is extremely unlikely within the standard
$\Lambda$CDM model (Menanteau et al. 2012; Jee et al. 2014).

The cluster exhibits a binary mass structure (Jee et al. 2014), with the two
sub-clusters denominated northwestern (NW) and southeastern (SE), respectively
(see, for example, Figure 1 of Ng et al. 2015), because of their positions.
Weak lensing (WL) analysis (Jee et al. 2014) shows a projected separation
between the mass centroids of the two components of approximately $d\sim
700\rm{\,kpc}$, with a mass ratio in the range $\sim 2:1$. Dynamical estimates
(Menanteau et al. 2012) suggest for the merging a high relative velocity for
the merging infall, in the range between $\sim
1,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$ and $\sim
2,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$, with the uncertainties in the
allowed range being due to projection effects.

A significant effect expected in high-velocity merging clusters is the
presence of large spatial offsets between the X-ray surface brightness peaks
and the centroid of the surface mass densities (Molnar et al. 2012), as well
as between the SZ and X-ray centroids. All of these features have been
observed in the El Gordo cluster (Menanteau et al. 2012; Zitrin et al. 2013):
the measured distance between the SZ and X-ray centroid is $d\sim
600\rm{\,kpc}$ while the SZ centroid is offset by a distance of $\sim
150\rm{\,kpc}$ from the main mass centroid. Additionally, the peak of the
X-ray emission is also offset from the position of the Brightest Cluster
Galaxy (BCG, Menanteau et al. 2012). This is located in the SE cluster
component, and with respect the SE X-ray peak is offset by $\sim
70\rm{\,kpc}$.

Moreover, the X-ray morphology of El Gordo presents several interesting
features. The total X-ray luminosity is $L_{X}\sim 2\cdot
10^{45}\rm{\,ergs^{-1}}$ in the $0.5-2$ keV band (Menanteau et al. 2012), with
the highest X-ray emission being in the SE region and being characterized by a
single peak and two extended faint tails elongated beyond the peak. The X-ray
emission in the NW region is significantly weaker and presents a wake-like
feature, indicative of turbulent flows. This well defined X-ray morphology
suggests that the merging axis is quite close to the plane of the sky
(Menanteau et al. 2012), this is also consistent with the observed presence of
a double radio relic (Lindnet et al. 2014).

This set of multifrequency observations leads to the the conclusion that El
Gordo is a cluster having a major merger at high redshift, which consists of
two massive clusters colliding at a high velocity ($\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}2,000\rm{\,km\leavevmode\nobreak\
s^{-1}}$). In this context the most likely physical interpretation of the
ongoing merging (Menanteau et al. 2012; Jee et al. 2014) is a scenario in
which the dense cool gas core of the secondary is moving from NW to SE,
inducing in the ICM of the primary a wake and the two-tail cometary structure
seen in the X-ray images. The originally present gas core of the primary is
then absent because it has been destroyed during the collision with the high-
velocity compact core of the secondary. The latter is then moving away from
the primary and the two clusters are seen after the first pericenter passage,
this is the so called outgoing scenario.

These observational features indicate that the El Gordo cluster can be
considered an extreme merging event like the Bullet cluster, but at a higher
redshift. High-velocity collisions between massive clusters constitute a
severe challenge to the commonly accepted $\Lambda$CDM paradigm, for example
Asencio et al. (2021) argue that in a $\Lambda$CDM context the detection of a
cluster like El Gordo can be ruled out with a high significance.

Another noteworthy feature of El Gordo is the peak location of the different
mass components, as probed by various multifrequency observations (Menanteau
et al. 2012; Jee et al. 2014; Zitrin et al. 2013; Diego et al. 2020; Kim et
al. 2021). As expected, the X-ray peak position of the SE cluster is spatially
offset from that of its mass centroid; however the X-ray peak is not trailing
the DM peak as in the Bullet Cluster but, on the contrary, is leading the mass
peak. This is clearly at variance with what is expected from dissipative
arguments, but it can be explained in the so called returning scenario, as
proposed by Ng et al. (2015). According to those authors, the two clusters are
not moving away from each other, i.e. just after the pericenter passage, but
the cluster is now in a post apocenter phase with the two clusters moving
towards each other.

Moreover the BCG is not only trailing the X-ray peak as well, but it also
appears to be spatially offset from the mass centroid. This behavior has been
observed in other clusters (Hallman et al. 2004; Canning et al. 2012), and in
an SIDM scenario it is expected to exist in major mergers, between the
collisional DM and collisionless galaxies (Kim et al. 2017).

Because of its many unusual properties, El Gordo has been the subject of
several N-body/hydrodynamical simulations aimed at reproducing the observed
X-ray morphology (Donnert 2014; Molnar & Broadhurst 2015; Zhang et al. 2015,
2018). The simulations of Donnert (2014) were able to reproduce some of the
cluster X-ray properties, such as the total X-ray luminosity and the offset
between the X-ray peak the SZ centroid, but failed to replicate the observed
twin-tailed X-ray morphology. This could be due to a number of reasons, but
the numerical scheme used in the simulations could be playing a role. Donnert
(2014) carried out the simulations by employing a smoothed particle
hydrodynamics (SPH, Price 2012) scheme, which in its standard formulation it
is well known to suppress the growth of instabilities and, in turn, turbulent
mixing (Agertz et al. 2007; Mitchell et al. 2009).

For this reason, in their simulations of the El Gordo cluster, Molnar &
Broadhurst (2015) chose to use the Eulerian adaptative mesh refinement (AMR)
code FLASH (Fryxell et al. 2000). By considering a range of initial merger
parameters they were able to reproduce the twin-tailed X-ray morphology, but
with a projection angle between the plane of the sky and the merging axis
which is higher ($\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}40^{\circ}$) than that suggested by
radio relic polarization arguments (Ng et al. 2015). Moreover the total
cluster X-ray luminosity and mean temperature extracted from the simulations
were lower than the measured values, whilst the spatial offsets between the
centroids were larger than those observed.

The most complete simulation study of El Gordo presented so far was computed
by Zhang et al. (2015, hereafter ZYL15). By considering a wide range of
initial merging parameter space, the authors constructed a large set of
N-body/hydrodynamical binary merging simulations aimed at reproducing the
observational aspects of El Gordo. They used both SPH and AMR based codes.

From their simulation ensemble the model which is best able to reproduce many
of the observed properties of El Gordo is a merging cluster (‘model B’) with
total mass $\sim 3\cdot 10^{15}\rm{\,M_{\odot}}$ and a high mass ratio ($\sim
3.6$). The initial merging configuration is that of an off-axis merger with
impact parameter $\sim 800\rm{\,kpc}$ and a very high initial relative
velocity ($\sim 2,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$). The best match
with El Gordo measured offsets and X-ray images, as well as temperature and
luminosity, is obtained $t\sim 0.14$ Gyr after passing the pericenter and with
a modest inclination angle ($\sim 30^{\circ}$) between the merging axis and
the sky plane.

Model B of ZYL15 can match most of the observed features of the El Gordo
cluster, but it is not entirely free of problems. Specifically a final single
X-ray peak is found if the gas core of the primary is destroyed during the
collision, but this is obtained only if the initial gas fraction of the
primary is low ($\sim 5\%$). This in turn produces X-ray emission in the outer
region of the merging cluster which is lower than observed. Finally, it worth
noting that the simulations of ZYL15 can reproduce the total X-ray luminosity
of El Gordo only for a primary cluster mass $\sim 2.5\cdot
10^{15}\rm{\,M_{\odot}}$. To date, the simulation set of ZYL15 constitutes the
most exhaustive simulation study of the El Gordo cluster. However the merging
scenario investigated by ZYL15 might need to be revisited, in the light of
recent works suggesting substantially smaller masses for the El Gordo cluster
than those previously measured.
By adopting a free-form lens model, Diego et al. (2020) estimate from their SL
study a cluster total mass of the order of $\sim 10^{15}\rm{\,M_{\odot}}$.
This value is about a factor of two lower then previous SL studies (Zitrin et
al. 2013). In a recent paper Kim et al. (2021) presented an improved WL
analysis aimed at estimating the total mass of El Gordo and that of its
components. The values they reported for the NW and SE clusters are $\sim
9.9\cdot 10^{14}\rm{\,M_{\odot}}$ and $\sim 6.5\cdot 10^{14}\rm{\,M_{\odot}}$,
respectively. These values are significantly lower ($\sim 60\%$) than those
dynamically inferred by Menanteau et al. (2012). They are also lower ($\sim
30\%$) than the WL values reported by Jee et al. (2014).
These updated mass estimates have been obtained independently from studies by
both SL (Diego et al. 2020) and WL (Kim et al. 2021), based on improved
lensing modeling; they both give consistently lower mass values than those
estimated in previous works. This clearly has the advantage of easing the
tension with the $\Lambda$CDM model (Kim et al. 2021), but see the recent
study of Asencio et al. (2023) for a different view. Nonetheless, it remains
to be tested whether this range of masses for the El Gordo cluster is
consistent with its observed X-ray morphology.
Therefore the aim of this paper is to perform a series of
N-body/hydrodynamical simulations, in a merging scenario incorporating the
recently revised El Gordo cluster masses. The simulations are performed for a
variety of merging initial conditions: we do this in order to find the final
merging configuration which can best reproduce various observations of El
Gordo. The initial condition set up mirrors that implemented by ZYL15; this
has been purposely chosen in order to ease the comparison of our results with
previous findings.
The structure of the paper is the following. Section 2 briefly describes the
hydrodynamic method we use, as well as the procedures employed to set up the
merging initial conditions, illustrating how to construct particle
realizations of cluster halos initially at equilibrium and the initial merging
kinematics. The results are presented in Section 3, with Section 3.1 showing
results from off-axis merger simulations aimed at reproducing the observed
twin-tailed X-ray morphology. In Section 3.2 we analyze the consistency of a
returning scenario according to the results extracted from head-on cluster
collisions, while Section 3.3 is dedicated to merger simulations incorporating
a BCG stellar component and to the impact of self-interacting DM on some
merging models. Section 4 summarizes our main conclusions. Throughout this
work we use a concordance $\Lambda$CDM cosmology, with $\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$ and Hubble constant $H_{0}=70\equiv 100h$ km s-1 Mpc-1.

##  2 SIMULATION SETUP

###  2.1 Simulations

Our code consists of an improved SPH numerical scheme for the hydro part,
coupled with a standard treecode (Hernquist & Katz 1989) to solve the gravity
problem. The Lagrangian SPH code employs an entropy conserving formulation,
whilst in the momentum equations SPH gradients are estimated by evaluating
integrals and performing a matrix inversion.

This tensor approach (García-Senz et al. 2012) has been tested in a variety of
hydrodynamical test cases (Valdarnini 2016), showing that it outperforms
standard SPH (Price 2012) by greatly reducing gradient errors and in terms of
code performances can be considered competitive with results obtained from AMR
codes (Valdarnini 2016). The code has previously been applied to the study of
turbulence in galaxy clusters (Valdarnini 2019) and to that of cool-cores in
merging clusters (Valdarnini & Sarazin 2021). In particular, we refer to this
last study for a comparison with results extracted from merger simulations
based on Eulerian AMR schemes and for a full description of the code. In the
following we will refer to the SPH scheme used here as integral SPH (ISPH).

For the sake of completeness we introduce the SPH gas density equation, to
which we will later refer for further generalization. In SPH there are $N$ gas
particles with mass $m_{i}$, position $\bf r_{i}$ and velocity $\bf v_{i}$
which represent the fluid.

The SPH gas density $\rho({\bf r})$ at the particle position $\bf r_{i}$

$\rho({\bf r_{i}})=\sum_{j}m_{j}W(|{\bf r}_{ij}|,h_{i})\leavevmode\nobreak\ ,$
(1)

gives the particle gas density $\rho(\bf r_{i})\equiv\rho_{i}$. The summation
should be over all the set of $N$ particles, but because the adopted kernel
$W(|{\bf r}_{ij}|,h_{i})$ has compact support it reduces to only the $N_{nn}$
neighboring particles $j$ for which $|{\bf r}_{i}-{\bf r}_{j}|\leq 2h_{i}$. We
use here the $M_{4}$ kernel (Price 2012). The smoothing length $h_{i}$ is
obtained by finding the root of

$\frac{4\pi(2h_{i})^{3}\rho_{i}}{3}=N_{nn}m_{i}.$ (2)

for which we set $N_{nn}=32\pm 3$.

###  2.2 Initial conditions

To study the merging cluster El Gordo we perform a series of
N-body/hydrodynamical ISPH simulations. The merging system consists of a
primary cluster with mass $M_{1}$ and a secondary with mass $M_{2}$. The
latter is specified by the given mass ratio $q=M_{1}/M_{2}$ .

The cluster mass is defined as the mass $M_{200}$ within the cluster radius
$r_{200}$. This is the radius $r_{\Delta}$ within which the total mass is

$M_{\Delta}=\frac{4\pi}{3}\Delta\rho_{c}(z)r_{\Delta}^{3}\leavevmode\nobreak\
,$ (3)

and the mean density is $\Delta$ times the cosmological critical density
$\rho_{c}(z)$. The cluster radii are calculated by setting $z=0.87$, the
redshift of the El Gordo cluster.

The initial condition setup that we use here for the two colliding clusters
follows the same setting adopted by ZYL15. This is purposely chosen with the
intent of validating our initial condition method by comparing results
extracted from similar merging runs, and then comparing different merger
simulations against a baseline merging model.

To construct the initial conditions for our merger simulations we first
perform a particle realization of two individual spherical halos, initially in
hydrostatic equilibrium. Each halo consists of multiple mass components: dark
matter, gas and eventually star particles.

####  2.2.1  Dark matter

For the DM density profile we assume an NFW profile

$\begin{array}[]{llll}\rho_{DM}(r)&=&\dfrac{\rho_{s}}{r/r_{s}(1+r/r_{s})^{2}}\,,&0\leq
r\leq r_{200}\,,\\\ \end{array}$ (4)

where $\rho_{s}$ and $r_{s}$ are the scaling parameters for density and
radius, respectively. The scale radius $r_{s}$ is defined by
$r_{s}=r_{200}/c_{200}$, where $c_{200}$ is the concentration parameter. We
set the value of $c_{200}$ using the $c-M$ relation of Duffy et al. (2008)

$c_{200}=5.71a^{0.47}\left(\frac{M_{200}}{2\cdot 10^{12}\rm{\,M_{\odot}{\it
h}^{-1}}}\right)^{0.084}\leavevmode\nobreak\ ,$ (5)

where $a=1/(1+z)$. The DM cluster density profile is then specified once the
mass $M_{200}$ and the concentration parameter $c_{200}$ are given.

At radii larger than $r_{200}$ the cumulative mass profile diverges as
$r\rightarrow\infty$. We therefore introduce for $r>r_{200}$ an exponential
cutoff up to a final radius $r_{max}=\nu r_{200}$:

$\begin{array}[]{lll}\rho_{DM}(r)&=&\rho_{DM}(r_{200})(r/r_{200})^{\delta}\exp{-\displaystyle{\left(\frac{r-r_{200}}{r_{decay}}\right)}}\,,\\\
&&\leavevmode\nobreak\ \leavevmode\nobreak\
r_{200}<r<r_{max}\leavevmode\nobreak\ ,\end{array}$ (6)

where $r_{decay}=\eta r_{200}$ is the truncation radius, and the parameter
$\delta$ is set by requiring the first derivative of the DM density profile to
be continuous at $r=r_{200}$. For the runs presented here, we set the
truncation parameters to the values $(\nu,\eta)=(2,0.2)$; this choice being
motivated by stability considerations (see Valdarnini & Sarazin 2021, for a
discussion about the optimum choice of these parameters).

To construct a particle realization of the DM density profile we first compute
the cumulative DM mass $M_{DM}(<r)$ within the radius $r$. We solve for the
particle radius $r$ by inverting $q(r)=M_{DM}(<r)/M_{DM}(<r_{rmax})=y$, where
$y$ is a uniform random number in the interval $[0,1]$.

Particle speeds are determined according to Kazantzidis et al. (2004) . For
spherically symmetric systems the density $\rho_{m}(r)$ of a given mass
component is given by

$\rho_{m}(r)=4\pi\int_{0}^{\Psi(r)}f_{m}(\mathcal{E})\sqrt{2[\Psi(r)-\mathcal{E}]}\leavevmode\nobreak\
d\mathcal{E}\leavevmode\nobreak\ ,$ (7)

where $\Psi(r)=-\Phi(r)$ is the relative gravitational potential of the system
, $\mathcal{E}=\Psi-v^{2}/2$ is the relative energy and $f_{m}(\mathcal{E})$
the corresponding distribution function. Here the subscript $m=DM,s$ denotes
the mass (DM or stars) component under consideration.

Through the Abell integral transform (Binney & Tremaine 1987) it is possible
to invert Equation (7) to obtain

$f_{m}(\mathcal{E})=\frac{1}{\sqrt{8}\pi^{2}}\left[\int_{0}^{\mathcal{E}}\frac{d^{2}\rho_{m}}{d\Psi^{2}}\frac{d\Psi}{\sqrt{(\mathcal{E}-\Psi)}}\right]\leavevmode\nobreak\
,$ (8)

where the second order derivative ${d^{2}\rho_{m}}/{d\Psi^{2}}$ can be
evaluated analytically from $\rho_{m}(r)$ and $M(<r)$.

We construct tables of $f_{m}(\mathcal{E})$ by evaluating numerically the
integral (8) over a range of energies. We then use an acceptance–rejection
method to obtain the particle speed $v=\sqrt{2[\Psi(r)-\mathcal{E}]}$ for a
particle at position $r$.

Finally, the directions of the particle position and velocity vectors are
assigned by randomly orienting unit vectors.

Table 1:  IDs and initial collision parameters of the off-axis merger
simulations of Sections 3.1 . a

Model |  $M^{(1)}_{200}$ $[\rm{\,M_{\odot}}]$ | $r_{200}\leavevmode\nobreak\ [\rm{\,Mpc}]$ | $q$ | $V\leavevmode\nobreak\ [\rm{\,km\leavevmode\nobreak\ s^{-1}}]$ | $P\leavevmode\nobreak\ [\rm{\,kpc}]$ | $(f_{g1},f_{g2})$ | gas profile
---|---|---|---|---|---|---|---
B_1 | $2.5\times 10^{15}$ | 2.02 | 3.6 | 2500 | 800 | (0.05,0.1) | Burkert
A_1 | $1.3\times 10^{15}$ | 1.63 | 2.0 | 3000 | 300 | (0.1,0.1) | Burkert
Bf | $1.6\times 10^{15}$ | 1.74 | 2.32 | 2500 | 600 | (0.1,0.1) |  $\beta$-model
Bg | $1.6\times 10^{15}$ | 1.74 | 2.32 | 2000 | 600 | (0.1,0.1) |  $\beta$-model
Bh | $1.6\times 10^{15}$ | 1.74 | 2.32 | 1500 | 600 | (0.1,0.1) |  $\beta$-model
Bk | $1\times 10^{15}$ | 1.5 | 1.54 | 1500 | 600 | (0.1,0.1) |  $\beta$-model
Bl | $1\times 10^{15}$ | 1.5 | 1.54 | 2000 | 600 | (0.1,0.1) |  $\beta$-model

Notes. a Columns from left to right: ID of the merging model, halo mass
$M^{(1)}_{200}$ of the primary , cluster radius $r_{200}$ at which
$\Delta=200$, primary-to-secondary mass ratio $q=M_{1}/M_{2}$, initial
collision velocity, collision impact parameter, primary and secondary cluster
gas mass fractions $f_{g}$ at $r_{200}$, adopted model for the gas density
profile of the primary. .

####  2.2.2  Baryonic matter

To construct the cluster gas initial conditions we assume spherical symmetry
and hydrostatic equilibrium. Following ZYL15, for the gas density we adopt the
Burkert (1995) density profile

$\begin{array}[]{llll}\rho_{gas}(r)&=&\dfrac{\rho_{0}}{(1+r/r_{c})\left[1+r/r_{c}\right]^{2}}\,,&0\leq
r\leq r_{200}\,,\\\ \end{array}$ (9)

where $r_{c}$ is the gas core radius and $\rho_{0}$ the central gas density.
At radii larger than $r_{200}$ the density profile is continued by assuming
that it follows that of the DM (Zhang et al. 2014):

$\rho_{gas}(r)=\rho_{DM}(r)\frac{\rho_{gas}(r_{200})}{\rho_{DM}(r_{200})},\leavevmode\nobreak\
\leavevmode\nobreak\ r_{200}<r<r_{max}\leavevmode\nobreak\ .$ (10)

The choice of the core radii is observationally motivated by the X-ray
morphology of El Gordo. In their paper, ZYL15 chose to set $r_{c}=r_{s}/2$ and
$r_{c}=r_{s}/3$ for the primary and secondary cluster, respectively. These
choices will be discussed later, here for the initial gas density profile of
the primary we additionally consider the commonly employed non-isothermal
$\beta$-model (Cavaliere & Fusco-Femiano 1978)

$\rho_{gas}(r)=\rho_{0}\left(1+\frac{r^{2}}{r_{c}^{2}}\right)^{-\frac{3}{2}\beta}\leavevmode\nobreak\
.$ (11)

For a given gas core radius $r_{c}$ and slope parameter $\beta$, to completely
specify the gas density profile it is necessary to determine the central
density $\rho_{0}$. This is done by assigning the initial gas mass fraction
$f_{g}$ of the cluster at $r=r_{200}$ and then finding the corresponding root
value $\rho_{0}$.

The choice of the initial gas density profiles is a crucial issue in shaping
the final X-ray morphology of the merging system, we therefore postpone to a
dedicated Section (2.5) a thorough discussion of the adopted profiles and
their motivation.

Finally, under the assumption of hydrostatic equilibrium the gas temperature
at the radius $r$ can now be derived (Mastropietro & Burkert 2008; Donnert
2014)

$T(r)=\frac{\mu
m_{p}}{k_{B}}\frac{G}{\rho_{gas}(r)}\int_{r}^{\infty}\frac{\rho_{gas}(u)}{u^{2}}M_{tot}(<u)du\leavevmode\nobreak\
,$ (12)

where $M_{tot}(<r)$ is the total mass within the radius $r$ and we assume
$\mu=0.59\leavevmode\nobreak\ $. The gas thermal energy in the SPH equations
is obtained by setting $\gamma=5/3$ for the ideal equation of state.

The initial gas particle positions are generated starting from a uniform glass
distribution of points. The desired density profile is then obtained by
applying to the particle coordinates a radial transformation, such that the
corresponding mass profile is now consistent with the initial gas mass
profile. The initial gas velocities are set to zero.

####  2.2.3 Stellar component

For some of our merging runs we also consider halos which consist of DM, gas
and an additional star matter component. This is supposed to mimick the BCG
mass distribution, which we locate initially at the center of the halo.

For the density profile of the stellar component we use (Merritt et al. 2006;
ZuHone et al. 2019)

$\begin{array}[]{lll}\rho_{\star}(r)&=&\rho_{e}\exp{\bigl{\\{}-d_{n_{s}}\displaystyle{\bigl{[}(r/r_{e})^{1/n_{s}}-1\bigr{]}}\bigr{\\}}}\,,\\\
d_{n_{s}}&\sim&3n_{s}-1/3+0.0079/n_{s},\leavevmode\nobreak\
\leavevmode\nobreak\ {\it for}\leavevmode\nobreak\ n_{s}\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}0.5,\end{array}$ (13)

where $r_{e}$ is the effective radius and $\rho_{e}$ the stellar density at
that radius, we set here $n_{s}=6$. The BCG mass is derived according to the
relation of Kravtsov et al. (2018, their Table 2)

$\log_{10}M_{\star,BCG}\sim 0.38[\log_{10}(M_{500})-14.5].\leavevmode\nobreak\
$ (14)

For $M_{500}\sim 10^{15}\rm{\,M_{\odot}}$, which is appropriate for the range
of halo masses considered here, the relation (14) gives $M_{\star,BCG}\sim
2.3\cdot 10^{12}\rm{\,M_{\odot}}$. We can now obtain the density $\rho_{e}$ by
solving for the ratio $M_{\star,BCG}/M_{500}$ within $r_{500}$, once $r_{e}$
is given.

If one adopts the relation of Shen et al. (2003, Equation 33), then for
$M_{\star,BCG}\sim 3\cdot 10^{12}\rm{\,M_{\odot}}$ one obtains $r_{e}\simeq
30$ kpc. The smaller is the value of the effective radius, the higher must be
the central density for a given $M_{\star,BCG}$. However the mass resolution
of the simulations sets a lower limit for the minimum spatial scale which can
be resolved, so that the value of $r_{e}$ should be considered simulation
dependent. For example in their simulations ZuHone et al. (2019) set
$r_{e}=175$ kpc, as a compromise we set here $r_{e}=60$ kpc. This value is a
bit higher than the adopted gravitational softening length for the DM
particles (see below). For a $\sim 10^{15}\rm{\,M_{\odot}}$ halo mass and
$M_{\star,BCG}\sim 3\cdot 10^{12}\rm{\,M_{\odot}}$, we then obtain for the
central value of the stellar mass density $\rho_{\star}(r\sim
10^{-2}r_{200})\simeq 10^{4}\rho_{c}\simeq 3.7\cdot
10^{6}\rm{\,M_{\odot}}kpc^{-3}$. For a given initial star density distribution
we then determine positions and velocities of the star particles according to
the procedure described in Section 2.2.1.

Finally, we assign the mass of DM and gas particles as in the merger
simulations of Valdarnini & Sarazin (2021):

$m_{d}\simeq 8\times 10^{8}\rm{\,M_{\odot}}(M^{(2)}_{200}/2\times
10^{14}\rm{\,M_{\odot}})\leavevmode\nobreak\ $ (15)

and $m_{g}=f_{b}m_{d}/(1-f_{b}))\simeq 0.16m_{d}$, respectively. Here
$M^{(2)}_{200}$ is the mass of the secondary and $f_{b}\simeq 0.16$ the
cosmological mass gas fraction. With this mass resolution we have
$N_{DM}\simeq 3.4\times 10^{5}$ DM particles for an halo mass of $M_{200}\sim
10^{15}\rm{\,M_{\odot}}$, and $N_{g}\simeq 1.7\times 10^{5}$ gas particles for
an halo gas mass of $M_{g}\sim 10^{14}\rm{\,M_{\odot}}$.

Table 2: IDs of the merger simulations with initial merging parameters as
given by the corresponding merger models of Table 1. For each merger model the
additional subscripts refer to simulations with an initially different gas
core radius $r_{c}$ for the primary. For each simulation we report the value
of $r_{c}$ in kpc and that of the ratio $\zeta=r_{s}/r_{c}$, where $r_{s}$ is
the NFW scale radius of the primary. The units of the collision parameters
$\\{M_{1},r_{s},\leavevmode\nobreak\ q,\leavevmode\nobreak\
P,\leavevmode\nobreak\ V\\}$ for each model are the same as for Table 1, the
NFW radius $r_{s}$ is given in Mpc.

Model: $\\{M_{1},r_{s},\leavevmode\nobreak\ q,\leavevmode\nobreak\ P,\leavevmode\nobreak\ V\\}$ | merger simulation ID
---|---
|  |  |  |  |  |
Bf: $\\{1.6\cdot 10^{15},0.696,\leavevmode\nobreak\ 2.32,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 2,500\\}$ | Bf$\\_$rc29 | Bf$\\_$rc26 | Bf$\\_$rc22 | Bf$\\_$rc20 | Bf$\\_$rc17 | Bf$\\_$rc14
$r_{c}$ ($\zeta=r_{s}/r_{c}$) | 290 (2.4) | 260 (2.67) | 217 (3.2) | 200 (3.44) | 174 (4.) | 145 (4.81)
Bg: $\\{1.6\cdot 10^{15},0.696,\leavevmode\nobreak\ 2.32,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 2,000\\}$ | Bg$\\_$rc29 | Bg$\\_$rc23 | Bg$\\_$rc20 | Bg$\\_$rc17 | Bg$\\_$rc14 |
$r_{c}$ ($\zeta=r_{s}/r_{c}$) | 290 (2.4) | 232 (3) | 200 (3.437) | 174 (4) | 145 (4.81) |
Bh: $\\{1.6\cdot 10^{15},0.696,\leavevmode\nobreak\ 2.32,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 1,500\\}$ | Bh$\\_$rc29 | Bh$\\_$rc22 | Bh$\\_$rc20 | Bh$\\_$rc17 | Bh$\\_$rc16 | Bh$\\_$rc14
$r_{c}$ ($\zeta=r_{s}/r_{c}$) | 290 (2.4) | 217(3.2) | 200(3.44) | 174(4.) | 160(4.37) | 145(4.81)
Bk: $\\{1.\cdot 10^{15},0.574,\leavevmode\nobreak\ 1.54,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 1,500\\}$ | Bk$\\_$rc19 | Bk$\\_$rc18 | Bk$\\_$rc17 | Bk$\\_$rc15 | Bk$\\_$rc14 | Bk$\\_$rc12
$r_{c}$ ($\zeta=r_{s}/r_{c}$) | 191 (3.) | 179 (3.206) | 167 (3.43) | 155 (3.7) | 143 (4) | 119 (4.81)
Bl: $\\{1.\cdot 10^{15},0.574,\leavevmode\nobreak\ 1.54,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 2,000\\}$ | Bl$\\_$rc24 | Bl$\\_$rc21 | Bl$\\_$rc19 | Bl$\\_$rc18 | Bl$\\_$rc17 | Bl$\\_$rc14
$r_{c}$ ($\zeta=r_{s}/r_{c}$) | 240 (2.39) | 215 (2.67) | 191 (3.0) | 179 (3.207) | 167 (3.43) | 143 (4)

Previous tests (Valdarnini 2019; Valdarnini & Sarazin 2021) showed that with
these mass assignments the simulation numerical resolution is adequate to
properly describe the development of hydrodynamic instabilities during the
merging phase.

The gravitational softening parameters of the particles are set according to
the scaling (Valdarnini & Sarazin 2021)

$\varepsilon_{i}=15.8\cdot(m_{i}/6.2\times
10^{8}\rm{\,M_{\odot}})^{1/3}\rm{\,kpc}.$ (16)

For those merger simulations with a BCG, we decided to increase the number of
star particles $N_{s}$ by a factor of $4$. To obtain the mass of star
particles we then apply Equation (15), but with a coefficient scaled down by a
factor $4$. This is to prevent during the simulations possible discreteness
effects on the star particle distribution. With this choice of parameters the
lower limit $\varepsilon_{\star}\sim r_{200}/\sqrt{N_{s}}$ (Power et al. 2003)
to the gravitational softening length of the star particles is then consistent
with the values given by Equation (16). Note however that this lower limit is
derived by requiring discreteness effects to be negligible over a Hubble time,
whereas here the time scale of the collisions is $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1$ Gyr.

Particle positions and velocities are integrated using a standard hierarchical
block timestep scheme (Hernquist & Katz 1989). According to several accuracy
and stability requirements (Springel 2005), each particle $i$ has its own
optimal timestep $\Delta t_{i}$ which is chosen from a power-of-two hierarchy
$\Delta t_{m}=\Delta t_{0}/2^{m}$, where $m=0,1\ldots$ and $\Delta
t_{0}=1/100$ Gyr is the largest allowed timestep. During the simulations we
record particle properties at simulation times $t_{n}=n\Delta t_{0}$, where
$n=0,1,\ldots$ is the simulation step.

####  2.2.4 Initial merger kinematics

To construct the orbits of our merger simulations, we choose a Cartesian
system of coordinates $\\{x,y,z\\}$, with the mergers taking place in the
$\\{x,y\\}$ plane and the center of mass of the two clusters being centered at
the origin.

As previously indicated, the initial separation and relative velocity vectors
$\\{{\bf X}^{in},{\bf V}^{in}\\}$ of the two clusters are set as in ZYL15. For
a given mass ratio $q=M_{1}/M_{2}$ the initial separations are then ${\bf
X}_{1}=\\{d_{ini}/(1+q),P/(1+q),0\\}$ and ${\bf
X}_{2}=\\{-d_{ini}q/(1+q),-Pq/(1+q),0\\}$ for the primary and secondary
cluster, respectively. Here $P$ is the impact parameter of the collision and
$d_{ini}=2(r^{1}_{200}+r^{2}_{200})$. Similarly, for the initial velocities we
set ${\bf V}_{1}=\\{-V/(1+q),0,0\\}$ for the primary and ${\bf
V}_{2}=\\{Vq/(1+q),0,0\\}$ for the secondary, where $V$ is the initial
collision velocity.

We then shift particle positions and velocities of the two halos according to
the vectors $\\{{\bf X}^{in},{\bf V}^{in}\\}$, the merger dynamical evolution
being fully determined by the merging parameters
$\\{M_{1},\leavevmode\nobreak\ q,\leavevmode\nobreak\ P,\leavevmode\nobreak\
V\\}$. In terms of these parameters the energy of the system can be expressed
as

$E=\frac{1}{2}M_{1}V^{2}\left[\frac{1}{1+q}-\frac{GM_{1}}{d_{ini}}\frac{1}{V^{2}}\frac{2}{q}\right]\leavevmode\nobreak\
,$ (17)

where for the merging system we have assumed the absence of bulk motion and
that the initial distance $d_{ini}$ between the two clusters is greater than
the sum of their radii $r_{200}$.

###  2.3 Numerical implementation of Self-Interacting Dark Matter

The possibility that DM has collisional properties has been investigated by
many authors and is observationally motivated by the well-known challenges
faced by the standard CDM model over a wide range of scales. We refer to Tulin
& Yu (2018) for a recent review on the subject. On cluster scales, self-
interacting DM offers a natural explanation for the apparent spatial offsets
seen in merging clusters between the galactic and the DM component, though it
must be emphasized that this effect is not present in every merging cluster.

Several approaches have been proposed to implement DM self-interactions in
N-body simulations. A DM simulation particle is a ‘super’ particle
representing a patch of DM phase-space which contains a large number of
physical DM particles (Hockney & Eastwood 1988). The DM self-interactions are
then modeled according to a standard Monte Carlo method, in which a scattering
between two DM simulation particles takes place whenever a local scattering
probability is smaller than a $[0,1]$ random number.

In the following we restrict ourselves to the simplest case of isotropic and
elastic scattering between DM particles, moreover we further assume a
constant, velocity-independent, DM cross section $\sigma_{DM}$. Although an
SIDM model with a constant DM cross section fails to satisfy a number of
constraints over a wide range of mass scales (Kaplinghat et al. 2016),
nonetheless it can be considered a useful first-order approximation as long as
SIDM simulations are restricted to the range of cluster scales (ZuHone et al.
2019).

To determine the local scattering probability, different methods have been
proposed. The approach which we follow here has been introduced by
Vogelsberger et al. (2012), and its implementation requires the definition of
a local DM density.

For each DM particle $i$ we define a DM density $\rho_{DM}({\bf r_{i}})$
according to the SPH formulation, and as in Equation (1) we thus define a DM
smoothing length $h_{i}^{DM}$ such that

$\rho_{DM}({\bf r_{i}})=\sum_{j}m^{DM}_{j}W(|{\bf
r}_{ij}|,h^{DM}_{i})\leavevmode\nobreak\ ,$ (18)

where $m^{DM}_{i}$ denotes the mass of the DM simulation particle $i$ and the
summation is over $N_{nn}=32\pm 3$ DM neighboring particles. The DM smoothing
length $h_{i}^{DM}$ is evaluated whenever the particle $i$ is active and its
timestep is synchronized with the current simulation time. For a DM particle
which is not active (passive), its smoothing length is extrapolated to the
current time from its last evaluation.

For each DM particle $i$ the local scattering probability with a neighboring
DM particle $j$, within the simulation timestep $\Delta t_{i}$ is then
determined as (Vogelsberger et al. 2012)

$P_{ij}=m^{DM}_{i}W(r_{ij},h^{DM}_{i})\frac{\sigma_{DM}}{m_{X}}v_{ij}\Delta
t_{i}\leavevmode\nobreak\ ,$ (19)

where $v_{ij}=|{\bf v_{i}}-\bf{v_{j}}|$ is the relative velocity between
particles $i$ and $j$, ${m_{X}}$ is the physical mass of the DM particle and
$\Delta t_{i}$ the timestep of particle $i$.

The total probability for particle $i$ to scatter is then given by
$P_{i}=\sum_{j}P_{ij}/2$, where the factor of two accounts for the other
member of the scattering pair. A collision of particle $i$ with one of its
neighbors will occur if $P_{i}\leq x$, where x is a uniform random number in
the range $[0-1]$. According to Vogelsberger et al. (2012) the scattering
neighbor is selected by sorting the set of neighboring particles according to
their distance $r_{ij}$ from particle $i$. The scattering neighbor is then the
first particle $k$ which satisfies $x\leq\sum_{j}^{k}P_{ij}$.

When the collision condition is satisfied, in the case of isotropic scattering
the post-scattering velocities of the tagged pair are

$\left\\{\begin{array}[]{lll}{\bf u_{i}}&=&{\bf V}+(v_{ij}/2){\bf e}\\\ {\bf
u_{j}}&=&{\bf V}-(v_{ij}/2){\bf e}\leavevmode\nobreak\ ,\end{array}\right.$
(20)

where ${\bf V}=({\bf v_{i}}+{\bf v_{j}})/2$ is the center-of-mass velocity,
and $\bf e$ is a unit vector randomly oriented. Note that in the simulations
we set all the DM particles as having the same mass.

The scattering procedure described here is generically valid for a serial
code, however its implementation in a parallel code is not so straightforward
to arrange and it is described in Appendix A. We also test the scattering
implementation by showing results from an isolated SIDM halo.

Table 3: Initial merging parameters of a set of additional off-axis merger simulations with smaller impact parameters. The meaning of the parameters is the same as in Tables 1 and 2. Model |  $M^{(1)}_{200}$ $[\rm{\,M_{\odot}}]$ | $q$ | $V[\rm{\,km\leavevmode\nobreak\ s^{-1}}]$ | $P[\rm{\,kpc}]$ |  $r_{c}$[kpc] ($\zeta$)
---|---|---|---|---|---
Bja$\\_$rc14 | $1\times 10^{15}$ | 1.54 | 1500 | 400 | 143 (4)
Bjb$\\_$rc14 | $1\times 10^{15}$ | 1.54 | 2000 | 400 | 143 (4)
Bjc$\\_$rc12 | $1\times 10^{15}$ | 1.54 | 1500 | 500 | 119 (4.81)

Finally, in all of the merger simulations in which DM is self-interacting the
DM halos are assumed to follow initially the same truncated NFW profiles of
the $\sigma_{DM}=0$ simulations. We justify this assumption by assuming to be
subdominant the impact of initially cored DM density profiles on the evolution
of post-pericenter gas structures (ZuHone et al. 2019).

###  2.4 Observables and projection analyses

For each of our merger simulations we save particle properties at various
output times. This is so as to generate for a given viewing direction and at
the chosen epoch images of the surface mass density, X-ray surface brightness
and SZ amplitude. To validate our merger models, these maps can then be
compared against analogous images presented in previous works. In particular,
to consistently compare the mock maps extracted from our simulations with
those shown by ZYL15, we adopt for the observer the same coordinate system
$\\{\hat{x},\hat{y},\hat{z}\\}$ introduced by the authors to construct their
projected maps.

In such a system the line-of-sight is along the $\hat{z}$ direction and the
$\\{\hat{x},\hat{y}\\}$ coordinates refer to the cluster merger as seen in the
plane of the sky. The observer frame is obtained from the simulation frame
$\\{x,y,z\\}$ by applying two rotation matrices. The first performs a rotation
around the $z$-axis by an angle $\alpha$, the second rotates the tilted frame
by an angle $i$ around the $x^{\prime}$-axis. The angle $i$ between the axes
$z$ and $\hat{z}$ is thus the angle between the merging and sky plane. Our
viewing direction is thus defined by the angles $\\{\alpha,i\\}$, we refer to
ZYL15 for more details on the setup of the rotation angles.

For a given output time and viewing direction we integrate along the line-of-
sight to extract from the simulation particles 2D maps of observational
quantities. Specifically, we define the surface mass density as

$\Sigma_{m}(x,y)=\int_{los}\left[\rho_{gas}({\bf x})+\rho_{DM}({\bf
x})\right]dz\leavevmode\nobreak\ ,$ (21)

where $\rho_{gas}({\bf x})$ and $\rho_{DM}({\bf x})$ are defined at the
position $\bf x$ according to Equations (1) and (18), respectively.

To properly compare the observed Chandra X-ray images (Menanteau et al. 2012)
with the mock X-ray maps extracted from their simulations, ZYL15 applied the
MARX software package 111http://space.mit.edu/ASC/MARX to their input X-ray
maps. Here we follow a simpler approach and , as in Molnar & Broadhurst
(2015), we evaluate the X-ray surface brightness by integrating the X-ray
emissivity $\varepsilon(\rho_{g},T_{g},Z,\nu)$ along the line-of-sight and
over the energy

$\begin{split}\Sigma_{X}(x,y)=&\frac{1}{4\pi(1+z)^{4}}\int_{los}\,dz\\\
&\int\,\varepsilon(\rho_{g},T_{g},Z,\nu)A_{eff}(\nu)\,d\nu\leavevmode\nobreak\
.\end{split}$ (22)

where $T_{g}$ is the gas temperature and $\nu$ the frequency, $Z$ is the metal
abundance of the gas and $A_{eff}(\nu)$ the effective area of the telescope.
For the latter we use the 4 chip averaged ACIS-I effective area (Zhao et al.
2004).

However, unlike in Molnar & Broadhurst (2015), here the X-ray emissivity is
not calculated analytically but using cooling tables taken from the MAPPINGS
III library (Allen et al. 2008). As in ZYL15 we set the gas metallicity to
$Z=0.3Z_{\odot}$, with $Z_{\odot}\sim 0.02$ (Anders & Grevesse 1989) being the
solar value.

Finally, to consistently compare our maps with observations and previous
findings (ZYL15), we set the exposure time to $t_{exp}=60ks$ for all the
simulated maps. Our X-ray maps (22) will then be expressed in counts arc sec-2
and evaluated in the energy range $[0.5-2]$ keV, without taking into account
the impact of the cosmic X-ray background.

We define X-ray temperatures according to the adopted weighting scheme:

$T_{W}=\frac{\int T({\bf x}){\it W}d^{3}x}{\int{\it
W}d^{3}x}\leavevmode\nobreak\ ,$ (23)

where $W$ is the weight function. We consider for $\it W({\bf x})$ the
commonly employed gas density function ${\it W}=\rho_{\mathrm{g}}$ (mass-
weighted temperatures, $T_{\mathrm{mw}}$) and the weight function ${\it
W}=\rho^{2}_{\mathrm{g}}T^{-3/4}$, which defines the corresponding
spectroscopic-like temperature $T_{\mathrm{X}}$. This choice of the weight
function has been shown (Mazzotta et al. 2004) to provide an accurate
approximation of spectroscopic temperatures extracted from X-ray observations.
As in previous integrals, projected X-ray temperature maps are obtained by
integrating weighted temperatures along the the line-of-sight.

The SZ surface brightness is calculated at the frequency $\nu$ including
relativistic corrections (Itoh et al. 1998):

$\begin{split}\Sigma_{SZ}(x,y)=&\frac{\sigma_{T}k_{B}}{m_{e}c^{2}}\int_{los}n_{e}T_{g}\\\
&\left[g(\nu)+\Sigma_{k=1}^{k=4}Y_{k}\Theta^{k}\right]dz\leavevmode\nobreak\
,\end{split}$ (24)

where $m_{e}$, $k_{B}$, $\sigma_{T}$, $c$ and $n_{e}$ are the electron mass,
the Boltzmann constant, the Thomson cross section, the speed of light and the
electron number density, respectively. The function
$g(\nu)=\coth(x_{\nu}/2)-4$ is the non-relativistic frequency function, where
$x_{\nu}=h_{P}\nu/(k_{B}T_{cmb})$ and $T_{cmb}$ is the temperature of the
cosmic microwave background. The second term in the square brackets is a sum
of relativistic corrections, with the coefficients $Y_{n}$ being given in Itoh
et al. (1998) and $\Theta\equiv k_{B}T_{g}/m_{e}c^{2}$. As in ZYL15 we set
$\nu=150$ GHz and smooth the SZ maps with a Gaussian kernel with width
$\sigma_{SZ}=270$ kpc.

We evaluate maps of projected quantities on a 2D mesh of $N_{g}^{2}=512^{2}$
grid points with coordinates $\\{x_{g},y_{g}\\}$. Because of the Lagrangian
nature of SPH, the value of a function at the grid point ${\bf x}_{g}$ is the
sum

$f({\bf x}_{g})=\sum_{i}f_{i}\frac{m_{i}}{\rho_{i}}W(|{\bf x}_{g}-{\bf
x}_{i}|,h_{i})\leavevmode\nobreak\ $ (25)

over the subset of SPH particles which satisfy $|{\bf x}_{g}-{\bf x}_{i}|\leq
2h_{i}$. To perform the integrals along the line-of-sight we exploit the
property of the $M_{n}$ splines (Price 2012), which approximate a Gaussian in
the limit $n\rightarrow\infty$. For the $M_{4}$ kernel used here the integrals
along the line-of-sight $z$ are then performed analytically and we calculate
projected quantities at the grid point ${\bf x}_{g}$ as

$f({\bf
x}_{g})=\sum_{i}f_{i}\frac{m_{i}}{\rho_{i}}W_{2D}(x^{2}_{gi},\sigma^{G}_{i})\leavevmode\nobreak\
,$ (26)

where $x^{2}_{gi}=(x_{g}-x_{i})^{2}+(y_{g}-y_{i})^{2}$, $\sigma_{i}^{G}\simeq
h_{i}/\sqrt{3}$ (Dehnen & Aly 2012) and $W_{2D}$ is the 2D Gaussian

$W_{2D}=\frac{1}{(\sqrt{2\pi}\sigma^{G}_{i})^{2}}exp[-x^{2}_{gi}/2\sigma_{i}^{G}]^{2}.$
(27)

To locate the centroid positions of the various maps that we produce to study
the El Gordo cluster, we apply to the simulation particles the shrinking
circle method. This a 2D version of the shrinking sphere method, which is
commonly employed in simulations to locate DM density peaks (Power et al.
2003) or X-ray temperature peaks (Valdarnini 2006).

We describe the method in three dimensions, but the generalization to
projected quantities is straightforward. We first introduce the concept of
weighted mean position for a subset of simulation particles. This subset can
consist of either gas, DM or star particles which were initially members of
one of the two colliding clusters. We define a weighted mean position by
calculating for the particle subset the averaged center ${\bf
x_{c}}=\sum_{i}w_{i}{\bf x}_{i}/\sum_{i}w_{i}$ , where the weights $w_{i}$ are
chosen according to the centroid under consideration.

To locate the centroid position we then start from an initially large radius
$R^{(0)}(\sim 5\rm{\,Mpc})$ which contains all of the particle subset and we
calculate the weighted mean position $\@vec{x}_{c}^{(0)}$.

We then proceed iteratively to calculate at each step $k=1,2,..,M$ the
centroid position $\@vec{x}_{c}^{(k)}$, based on those subset particles which
are within a sphere of radius $R^{(k)}=fR^{(k-1)}<R^{(k-1)}$, located at
$\@vec{x}_{c}^{(k-1)}$. The last iteration $M$ is such that there must be at
least a number $N_{M}$ of particles within $R^{(M)}$. The cluster center is
then defined as $\@vec{x}_{c}^{(M)}$. We found that the results are robust for
$N_{M}\sim 100$, but the choice of the shrinking factor $f$ is critical: it
must be chosen sufficiently large ($f\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}0.8$) to avoid that the algorithm
fails when there are multiple peaks.

The method is easily applied in 2D to the projected particle positions along
the viewing direction, so that for a given mock image we calculate the
position of the mass centroids, as well as the location of the X-ray emission
and SZ peaks.

###  2.5 Merger models

As previously explained in the Introduction, the aim of this paper is to
determine whether it is possible to construct merger models for the El Gordo
cluster which can consistently reproduce the observed X-ray morphology, as
well as many of its physical properties. However for these models a
fundamental constraint must be that of having cluster masses in line with
recent WL studies, and not as high as those employed in previous simulations
(Donnert 2014; Molnar & Broadhurst 2015; Zhang et al. 2015).

Our merger models will be validated by contrasting our simulation results
against the corresponding ones presented by ZYL15 for their fiducial model B.
As already outlined, this is the model which is able to match most of the
observed El Gordo physical properties. For this reason we adopt here the same
initial merging configuration, and extract X-ray maps from our merger
simulations when the projected distance between the mass centroids is
$d_{DM}\sim 700\rm{\,kpc}$ after the pericenter passage. In all of the
considered models the orientation angles of the merging plane are taken to be
the same as in model B of ZYL15:
$\\{\alpha,i\\}=\\{-90^{\circ},30^{\circ}\\}$.

To significantly simplify the collision parameter space of our merger models
we assume for the secondary the same Burkert gas density profile (9) as in
model B of ZYL15, with the gas core radius set to $r_{c}=r_{s}/3$. We justify
this assumption by noting from Figure 14 of Menanteau et al. (2012) that the
map of the deprojected electron density exhibits a peak value of $n_{e}\simeq
5\cdot 10^{-2}\rm{\,cm^{-3}}$ for the SE component. This value is very close
to the core value found for the electron number density of the secondary in
model B , at the post-pericenter epoch which fits the observed projected
distances between the mass centroids. We use then this finding to argue that
the gas structure of the SE component is adequately described here by adopting
the initial gas density profile of the secondary.

We then construct our ensemble of off-axis mergers by considering different
values of the primary mass $M^{(1)}_{200}$, but modifying the mass ratio
$q=M_{1}/M_{2}$ such that the mass of the secondary stays close to the
estimated value $M^{(2)}_{200}\sim 6.5\cdot 10^{14}\rm{\,M_{\odot}}$. In
accord with recent mass estimates (Diego et al. 2020; Kim et al. 2021) we
examine mergers with the following masses for the primary:
$M^{(1)}_{200}=1.6\cdot 10^{15}\rm{\,M_{\odot}}$ and
$M^{(1)}_{200}=10^{15}\rm{\,M_{\odot}}$.

For a chosen value of the primary mass we consider mergers whose initial
conditions differ in the choice of the initial collision velocity $V$ and
impact parameter $P$. Table 1 lists the different off-axis merger models which
we analyze, together with the corresponding merging parameters
$\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
P,\leavevmode\nobreak\ V\\}$; models A_1 and B_1 are the fiducial models A and
B of ZYL15. We will also be considering head-on mergers (models HO), but we
defer until Section 3.2 our description of the initial conditions adopted for
those models.

The adopted strategy for setting up the initial conditions of our merging runs
is based on the following result, found by ZYL15. When the authors considered
the possibility of a merger with initial conditions similar to those of model
B, but with a smaller primary mass ($M^{(1)}_{200}=1.6\cdot
10^{15}\rm{\,M_{\odot}}$, panel f in their Figure 5), they found a total X-ray
luminosity much smaller than that given by model B. Then they use this finding
to conclude that a massive merger ($M^{(1)}_{200}\simeq 2.5\cdot
10^{15}\rm{\,M_{\odot}}$) is required in order to obtain an X-ray luminosity
as high as that observed in the El Gordo cluster.

Here, we argue that a viable solution for obtaining consistency between the
cluster masses recently measured and observed luminosities is to consider the
possibility of having an initial gas density profile for the primary different
from the Burkert one (9) previously adopted.

Specifically, we employ the well-known $\beta$-model (11) to describe the
initial radial gas density profile of the primary. For a given central density
$\rho_{0}$, the radial behavior of the gas density can be modified by varying
the gas scale radius $r_{c}$ and the slope parameter $\beta$, with the aim of
constructing the initial conditions for a merging system which can satisfy the
observational constraints previously described.

We further simplify the parameter space of our merger models by setting
$\beta=2/3$, as already adopted in merger simulations for the Bullet Cluster
(Mastropietro & Burkert 2008) and, in particular, by Donnert (2014) in his
simulation study of the El Gordo cluster. For a given set of merging
parameters $\\{M^{(1)}_{200},\leavevmode\nobreak\ \xi,\leavevmode\nobreak\
P,\leavevmode\nobreak\ V\\}$, this choice leaves us with the gas core radius
$r_{c}$ of the primary being the only parameter needing to be set up for the
initial conditions of our merger models.

We then construct our ensemble of merger simulations by creating, for each of
the off-axis merger models listed in Table 1, a subsample of merging clusters
with different initial conditions. For a specified set of collision parameters
$\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
V,\leavevmode\nobreak\ P\\}$ we setup the initial conditions of our models by
considering a range of values for the core radius $r_{c}$. Table 2 lists the
merger models which we consider, the ID of each model being given by the ID of
the corresponding model in Table 1, with a subscript referring to the chosen
value of $r_{c}$. This is given in terms of the dimensionless parameter
$\zeta=r_{s}/r_{c}$. Hereafter, the initial gas core radius of the primary
will be denoted as $r_{c}^{\rm{prm}}$.

Finally, we also consider an additional list of merger models. The collision
parameters of these models are given in Table 3.

Figure 1:  X-ray surface brightness images extracted from simulations of
models B_1 and A_1, for different viewing directions. The box size is
$1.6\rm{\,Mpc}$ and time is in Gyr, with $t=0$ at the pericenter passage. In
each panel the epoch $t$ is when the projected distance $d_{DM}$ between the
mass centroids is approximately $d_{DM}\sim 700\rm{\,kpc}$, the orientation
angles $\\{\alpha,i\\}$ of the merging plane are defined in Section 2.4. In
each panel are overlaid the log-spaced contour levels of the projected X-ray
surface brightness (red) and mass density (white). From the inside to outside
the contour levels of the X-ray surface brightness and of the surface mass
density are: $(6.6,4.4,2.9,1.9,1.2)\cdot
10^{-1}\rm{\,counts\leavevmode\nobreak\ arcsec^{-2}}$ and $(5.6,3.1,1.8)\cdot
10^{-1}\rm{\,gr{\,cm^{-2}}}$, respectively. The crosses indicate the projected
spatial locations of the mass (green) and X-ray surface brightness (red)
centroids, the yellow cross showing the position of the SZ centroid. The maps
can be directly compared with the corresponding ones in Figure 1 and 3 of
ZYL15, the orientation along the viewing direction of the merging plane in the
sky being the same. In particular, the mid panels of models A_1 and B_1
correspond to the fiducial models A and B of ZYL15 see their Table 2. For
these models the X-ray luminosity in the $0.5-2$ keV band is given in units of
$10^{45}\rm{\,ergs^{-1}}$. For the fiducial model B_1 (top mid panel) the
dashed orange line shows the location and orientation of the spatial region
used to extract the X-ray surface brightness profile. This model is chosen as
the baseline model of the simulations.

##  3 Results

This Section is dedicated to the presentation of our main results, which we
divide into three subsections. In Section 3.1 we present results obtained by
performing off-axis merging simulations of the El Gordo cluster, but with a
mass of the primary in line with recent lensing estimates and lower than that
considered in previous works. In section 3.2 we investigate the possibility of
a returning scenario to explain the observed X-ray morphology and spatial
offsets between different centroids. Finally, in Section 3.3 some of the
merging runs of Section 3.1 are revisited by performing the simulations with a
BCG stellar component added to the initial halo configurations, and we also
consider the possibility of merging simulations within an SIDM framework.

###  3.1 Search for the optimal merger models

####  3.1.1 Code validation and physical processes leading to the formation of
the observed X-ray morphology

We first discuss simulation results obtained from the merging models A_1 and
B_1 of Table 1, the initial conditions of these models being the same as those
for the corresponding models of ZYL15. The purpose of these simulations is to
provide a sanity check of the initial condition setup adopted here for our
merging simulations, as well as for the analysis procedures presented in
Section 2.4. The hydrodynamic performances of the ISPH code have already been
tested in previous papers (Valdarnini 2016; Valdarnini & Sarazin 2021) against
results extracted from Eulerian-based AMR codes.

For these two models in Figure 1 we show, at different epochs, the projected
surface mass density and X-ray surface brightness maps. The bottom panels are
for model A_1 and top panels for model B_1, respectively. Each panel is for a
different viewing direction $\\{\alpha,i\\}$ of the merging plane, as in
ZYL15. The bottom panels (model A_1) can then be compared directly with those
in Figure 1 of ZYL15, while the top panels (model B_1) correspond to those of
their Figure 3. Each panel refers to the post-pericenter epoch when the
projected distance between the mass centroids is $d_{DM}\sim 700\rm{\,kpc}$.

To ease comparison with the ZYL15 findings, in each panel we show contour
levels of DM (white) and X-ray surface brightness (red) maps. Additionally, we
also show the spatial location of the DM (green), X-ray (red) and SZ (yellow)
centroids.

A visual inspection reveals that our test simulations are in substantial
agreement with the corresponding ones presented in ZYL15, with the mock X-ray
images of Figures 1 which exhibit the same morphological features seen in
Figures 1 and 3 of ZYL15. In particular, for model A_1 the bottom left
($\\{\alpha,i\\}=\\{-50^{\circ},0^{\circ}\\}$) and right
($\\{\alpha,i\\}=\\{0^{\circ},0^{\circ}\\}$) panels show the same X-ray
structures seen in Figures 1a and 1c of ZYL15: a strongly asymmetric tail
(Figure 1a) and a small asymmetric wake-like feature (Figure 1c). The bottom
mid panel ($\\{\alpha,i\\}=\\{-50^{\circ},75^{\circ}\\}$) of Figure 1
corresponds to the fiducial model A of ZYL15 (Figure 1b). The X-ray morphology
is more symmetric here and very similar to that seen in Figure 1b, with a
bullet-like shape.

Similarly, the top panels of Figure 1 can be compared with those displayed in
Figure 3 of ZYL15 for their model B. The top left
($\\{\alpha,i\\}=\\{-90^{\circ},0^{\circ}\\}$) and right
($\\{\alpha,i\\}=\\{-90^{\circ},60^{\circ}\\}$) panels show the existence of a
wake structure trailing the secondary, which is asymmetric and oriented either
to the left or to the right, respectively. This is the same behavior as seen
in Figure 3a and 3c of ZYL15. The top mid panel of Figure 1 corresponds to
Figure 3b of ZYL15 and to their fiducial model B. A comparison with this model
is particularly significant, since this is the model which best matches the
observations. In this panel ($\\{\alpha,i\\}=\\{-90^{\circ},30^{\circ}\\}$)
the X-ray morphology is clearly characterized by a well defined twin-tailed
structure, as seen in observations and in their Figure 3b.

The consistency between our model B_1, with viewing angles
$\\{\alpha,i\\}=\\{-90^{\circ},30^{\circ}\\}$, and the fiducial model B of
ZYL15 can be put in a more quantitative way by comparing the respective X-ray
luminosities $L_{X}$ and spectral temperatures $T_{\mathrm{X}}$. The total
X-ray luminosity in the $0.5-2$ keV band is evaluated by means of a standard
SPH estimator, using the cooling tables implemented in our code. For the
spectroscopic-like temperature $T_{\mathrm{X}}$ we use Equation (23).

Figure 2: As in Figure 1 but for merger models Bf and Bg. Each panel refers to
a simulation with a different gas core radius $r_{c}=r_{s}/\zeta$ of the
primary (see Table 2). In all of the panels maps are extracted along the same
viewing direction ($\\{\alpha,i\\}=\\{-90^{\circ},30^{\circ}\\}$) of the
baseline model B_1 in Figure 1. The thresholds and the spacing of the contour
levels are the same as those in Figure 1. The left panels are the chosen
fiducial models Bf_rc20 and Bg_rc20 for these merging runs.

Additionally, a bow shock will form in front of the secondary as it reaches
the pericenter and climbs its way toward the apocenter, this is because of the
lower ICM sound speed due to the decrease in the gas temperature. The ICM gas
will thus be compressed and shock-heated to higher temperatures, forming a
leading edge oriented perpendicularly to the secondary’s direction of motion
(Ricker 1998; Mastropietro & Burkert 2008; Machado & Lima Neto 2013; Molnar &
Broadhurst 2015; Sheardown et al. 2019; Moura et al. 2021; Chadayammuri et al.
2022).

The values that we obtain for these quantities are given in Table 4 and are
$L_{X}\simeq 2.86\cdot 10^{45}\rm{\,ergs^{-1}}$ and $T_{\mathrm{X}}\simeq
12\rm{\,keV}$, respectively. These values can be directly compared with the
corresponding ones reported in Table 2 of ZYL15: $L_{X}\simeq 2.05\cdot
10^{45}\rm{\,ergs^{-1}}$ and $T_{\mathrm{X}}\simeq 15\rm{\,keV}$. From the
same table one also has $\delta T\simeq-1.1\cdot 10^{-3}\kern
2.0pt{}^{\circ}{\kern-2.0ptK}$ for the central SZ decrement, while from Table
4 $\delta T\simeq-9\cdot 10^{-4}\kern 2.0pt{}^{\circ}{\kern-2.0ptK}$. The
largest relative difference is between the X-ray luminosities, which is
naively expected given the proportionality of $L_{X}$ to the squared density
and the associated uncertainties.

Figure 3: As in Figure 2 but for merger models Bh.

Given that these results have been obtained using two numerically different
hydrodynamic schemes, and the uncertainties associated with the use of
different estimators, we conclude that there is a fair agreement between the
two simulations.

We finally conclude that our merger model B_1 is in accord with the
corresponding fiducial model B of ZYL15. We therefore adopt model B_1, which
we will simply refer to as model B from here on, as the baseline model against
which to validate the simulation results obtained from our merger models.

Before we discuss our findings from the merger models we have previously
introduced, it is useful to look first at the processes leading to the
formation of the twin-tailed morphology seen in the X-ray structure of the
merger model B of Figure 1. The main collision parameters which determine the
X-ray structure of the colliding clusters are the relative collision velocity
$V$ and the collision parameter $P$, while keeping fixed the mass ratio and
the viewing direction.

For a given set of collision parameters $\\{V,P\\}$ the merger will be
characterized by an infall phase in which the secondary falls in the potential
well of the primary, with a growing ram-pressure as the secondary approaches
the pericenter and the relative velocity increases. This pressure will strip
material from the secondary, leading to the formation of a downstream tail.
Because the collision is off-center, during its motion through the ICM of the
primary, the secondary will begin to deviate from its original trajectory and
the tail will no longer be aligned with the direction of motion.

The gas morphology which is produced after the collision thus depends on the
adopted set of merging parameters $\\{M_{1},\leavevmode\nobreak\
q,\leavevmode\nobreak\ P,\leavevmode\nobreak\ V\\}$. If $P$ gets higher while
the other collision parameters are kept fixed, there will be a significant
weakening of the ram pressure during the collision. This in turn will produce
a strong asymmetric one-tailed gas structure, which is not observed; moreover
the primary would not have been destroyed during the collision. Similarly, if
$P$ approaches zero, the collision will become axisymmetric around the
collision axis, in contrast to the observed twin-tailed X-ray morphology (see
next Section). The observed asymmetric two-tail X-ray structure is therefore
expected to be produced for just a certain range of values of the impact
parameter $P$.

Lowering the infall velocity $V$ without modifying the other collision
parameters decreases the amount of gas removed by ram pressure stripping, and
the final result will be a single tail trailing the secondary. Similarly, if
$V$ gets higher the interaction time between the two clusters will be strongly
reduced and the resulting X-ray map will be strongly asymmetric (cf. Figure 5a
and 5b of ZYL15).

We thus deduce that for a given set of merging parameters
$\\{M_{1},\leavevmode\nobreak\ q,\leavevmode\nobreak\ P,\leavevmode\nobreak\
V\\}$ the post-collision gas structures are expected to be roughly reproduced
if one reduces proportionally both $P$ and $V$. The same reasoning applies if
the mass of the primary is decreased. Finally, projection effects can hide the
real X-ray morphology and produce overlapping structures along the line-of-
sight.

####  3.1.2 Off-axis merger models with mass of the primary
$M^{(1)}_{200}=1.6\cdot 10^{15}\rm{\,M_{\odot}}$: models Bf, Bg and Bh

We now use the above findings to investigate different merger models aimed at
reproducing the observed El Gordo X-ray morphology in a revisited mass
estimate framework. As can be seen from Table 1, we will consider two families
of off-axis merger models. The first family has the mass of the primary set to
$M^{(1)}_{200}=1.6\cdot 10^{15}\rm{\,M_{\odot}}$ (models Bf, Bg and Bh), while
$M^{(1)}_{200}=10^{15}\rm{\,M_{\odot}}$ for the second family (models Bk, Bl
and Bj). This set of models will be presented in Section 3.1.3. All of the
models have the gas fraction of the primary set to $f_{g1}=0.1$. In our
merging runs the highest value that we use for the mass of the primary is
$M^{(1)}_{200}=1.6\cdot 10^{15}\rm{\,M_{\odot}}$. This is about $\sim 50\%$
higher than the mass estimate ($\sim 9.9\cdot 10^{14}\rm{\,M_{\odot}}$)
reported by Kim et al. (2021) for the NW cluster, but still within
observational uncertainties.

We present our simulation results following the top down ordering of Table 1,
with merging initial conditions which progressively deviate from those of the
reference model B. The first model which we consider is model Bf, with
$P=600\rm{\,kpc}$ and $V=2,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$, followed
by models Bg ($V=2,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$) and Bh
($V=1,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$).

Model Bf is nothing else than the merger model presented by ZYL15 in their
Figure 5f. The authors rule out this model because it is unable to match the
required total X-ray luminosity. Our goal here is to assess the extent to
which a given merger model, with a specific set of merging initial conditions
$\\{M_{1},\leavevmode\nobreak\ q,\leavevmode\nobreak\ P,\leavevmode\nobreak\
V\\}$, can reproduce the observational features of the El Gordo cluster.

As introduced in Section 2.5, we adopt here the simulation strategy of
performing a set of merging simulations with initial conditions which differ
in the choice of the primary’s gas scale radius $r_{c}^{\rm{prm}}$, while
keeping fixed the other collision parameters of the chosen merger model. Table
2 reports the set of merging simulations specific to a given model, together
with the corresponding gas core radii.

The top panels of Figure 2 show, from the left to the right, the mock X-ray
images extracted from the runs of models Bf_rc17, Bf_rc29 and Bf_rc20,
respectively. The current epoch is defined when the projected distance
$d_{DM}$ between the mass centroids is approximately $d_{DM}\sim
700\rm{\,kpc}$, and as for the fiducial model B of Figure 1 the viewing angles
are $\\{\alpha,i\\}=\\{-90^{\circ},30^{\circ}\\}$. Hereafter, the same
projection direction will be used to extract X-ray surface brightness maps
from the off-axis merging simulations presented here.

For the Bf models, the NFW scale radius $r_{s}$ of the primary takes the value
of $r_{s}\simeq 0.7\rm{\,Mpc}$, so that the $r_{c}^{\rm{prm}}$ ranges from
$r_{c}^{\rm{prm}}\sim 145\rm{\,kpc}$ ($\zeta=4.81$: Bf_rc14) up to
$r_{c}^{\rm{prm}}\sim 290\rm{\,kpc}$ ($\zeta=2.4$: Bf_rc29). From the maps of
Figure 2 it appears then that a change of $r_{c}^{\rm{prm}}$, within a certain
range of values and under certain conditions for the primary’s gas density
profile, has the same impact on the final X-ray morphology of the merger as if
one were modifying the initial impact parameter $P$.

Figure 4: As in Figure 2 but for the merger models Bk and Bl. Model Bl_rc24 is
the chosen fiducial model for the Bl runs.

Specifically, model Bf_rc29 clearly illustrates the effect on the post-
collision X-ray structure of the merging when the gas core radius is
increased. At the observer epoch, the X-ray morphology clearly exhibits only
one tail trailing the secondary (top mid panel of Figure 2). In accordance
with the previous discussion about the dependence of the final X-ray
morphology on the collision parameters, we argue that for model Bf_rc29 the
final one-tail X-ray feature is a consequence of the reduced ram pressure
experienced by the secondary during the collision. This is in turn due to the
primary’s lower gas density because of the larger core radius.

If instead $r_{c}^{\rm{prm}}$ approaches small values, the opposite must hold.
The amount of gas stripped by the secondary during its motion through the
primary’s ICM will be much higher because of the increased ram pressure. As
can be seen from model Bf_rc17 (top left panel of Figure 2) the X-ray
morphology at the observer epoch now exhibits a strongly asymmetric two-tailed
structure, with one tail previously absent and now with a much stronger
emission than the other.

These findings suggest that for the collision parameters of the Bf family
models, the characteristic twin-tailed X-ray structure seen in the El Gordo
cluster can be reproduced with an appropriate choice of $r_{c}^{\rm{prm}}$.
Model Bf_rc20 (top right panel of Figure 2) is the model whose X-ray image is
best able to match the requested morphology: for this model
$r_{c}^{\rm{prm}}\sim 200\rm{\,kpc}$ and its twin-tailed morphology well
reproduces that of fiducial model B in Figure 1. Hereafter, we will denote
model Bf_rc20 as the fiducial model for the Bf family of merger models.

The same reasoning can be applied to the Bg models, for which the bottom
panels of Figure 2 show the mock X-ray maps extracted from three runs with
different values of the gas core radius of the primary, $r_{c}^{\rm{prm}}$.
For these models the initial collision velocity is lower
($V=2,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$) than in the Bf models, and
model Bg_rc20 ( $r_{c}^{\rm{prm}}\sim 200\rm{\,kpc})$ is now the fiducial
model which produces the optimal X-ray morphology. Note that despite the
initial velocity being lower than in model Bf_rc20, the gas core radius is the
same, but the X-ray morphology is a bit more asymmetric.

Figure 5: As in Figure 4 but for merger models Bj of Table 3. The meaning of
the symbols is the same, as well as the viewing direction.

Figure 3 shows for the Bh models X-ray maps extracted from three merger
simulations. For these models $V=1,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$
and it is interesting to note that all of the considered models, from
$r_{c}^{\rm{prm}}\sim 290\rm{\,kpc}$ ($\zeta=2.4$: Bh_rc29) down to
$r_{c}^{\rm{prm}}\sim 146\rm{\,kpc}$ ($\zeta=4.81$: Bh_rc14), produce strongly
asymmetric X-ray structures which are in contrast with the observed X-ray
morphology. This is clearly a consequence of the now much weaker ram pressure
due to the relatively low initial velocity between the two clusters.

We use these findings to argue that merger models of the El Gordo cluster,
with a primary mass of approximately $M^{(1)}_{200}\simeq 1.6\cdot
10^{15}\rm{\,M_{\odot}}$, require collision velocities at least as high as
$V\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}2,000\rm{\,km\leavevmode\nobreak\
s^{-1}}$ in order to be able to reproduce the observed X-ray morphology.

####  3.1.3 Off-axis merger models with mass of the primary
$M^{(1)}_{200}=10^{15}\rm{\,M_{\odot}}$: models Bk, Bl and Bj

For the family of merger models with mass of the primary
$M^{(1)}_{200}=10^{15}\rm{\,M_{\odot}}$, we have model Bk ($P=600\rm{\,kpc}$
and $V=1,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$) and Bl
($V=2,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$), plus the additional models Bj
of Table 3.

Figures 4 and 5 show the X-ray images constructed from simulation results of
the merger models Bk, Bl and Bj. The top panels of Figure 4 demonstrate that
the observed twin-tailed morphology cannot be reproduced by merger models with
initial collision velocity as low as $V=1,500\rm{\,km\leavevmode\nobreak\
s^{-1}}$ (Bk runs), regardless of the chosen value of the primary’s gas core
radius. This behavior is similar to that of model Bh (Figure 3), but now the
mass of the primary is lower.

For Bl models, the collision velocity is increased to
$V=2,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$. The bottom panels of Figure 4
show that this class of merger models can produce an X-ray image which
compares well visually with those previously used to define the list of
fiducial models. Model Bl_rc24 (bottom left panel of Figure 4) is the merger
model which is best able to match the observed X-ray morphology: for this
model $r_{s}\simeq 0.58\rm{\,Mpc}$ and $r_{c}^{\rm{prm}}\sim 240\rm{\,kpc}$ is
the corresponding value of its primary’s gas core radius. This model is added
to our list of fiducial models.

So far we have considered only merger models with initial impact parameter
$P=600\rm{\,kpc}$, apart from the fiducial model B being an exception
($P=800\rm{\,kpc}$) because of the high value of its primary’s mass. In order
to assess the viability of merger models with a smaller value of the impact
parameter we ran an additional set of models (Bj), with $P$ now ranging from
$P=500\rm{\,kpc}$ down to $P=400\rm{\,kpc}$. For each set of collision
parameters $\\{V,\leavevmode\nobreak\ P\\}$ we performed merging simulations
by considering different values of $r_{c}^{\rm{prm}}$. We show here results
extracted from the most significant merging runs, Table 3 lists the collision
parameters of these models. For all of them the mass of the primary is set to
$M^{(1)}_{200}=10^{15}\rm{\,M_{\odot}}$.

Figure 5 shows the X-ray maps extracted from the merger models Bj of Table 3.
From these maps it is possible to deduce that mergers with initial impact
parameter as low as $P=400\rm{\,kpc}$ are unable to generate the requested
X-ray morphology. The X-ray structure is strongly asymmetric, with only one
tail present if $V=2,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$ (Bjb_rc14 run)
and approaches that of a head-on collision when
$V=1,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$ (Bja_rc14 run). The asymmetric
one-tail structure is even more pronounced if one considers model Bjc_rc12 for
which $P=500\rm{\,kpc}$.

####  3.1.4 Comparison of the results with previous findings and observations.

In order to better constrain the fiducial merger models which have been
singled out by our set of simulations, for each of these models we now compare
results and observations. Table 4 lists the measured values of several
quantities, together with their observational estimates. These are taken from
Table 2 of ZYL15.

The spectral X-ray temperatures $T_{\mathrm{X}}$ given in Table 4 are obtained
by using the estimator given by Equation (23). This approach is much simpler
than that adopted by ZYL15, who estimated spectral temperatures by fitting
mock X-ray spectra using the software package MARX, nonetheless it has been
employed by many authors (Mastropietro & Burkert 2008; Molnar et al. 2013;
Zhang et al. 2018) and it provides a good approximation to cluster spectral
fit temperatures, as demonstrated by Nagai et al. (2007).

The spectral temperature of model B ($T_{\mathrm{X}}\simeq 12\rm{\,keV}$) is
smaller than that reported by ZYL15 for their fiducial model B
($T_{\mathrm{X}}\simeq 15\rm{\,keV}$) by about a factor of $\sim 20\%$. As
previously discussed we do not consider this discrepancy particularly
significant, given the use of different codes and estimators. The temperatures
of models Bf_rc20 and Bg_rc20 ($T_{\mathrm{X}}\simeq 13.5\rm{\,keV}$) fare
better than model B to approximate the observational estimate
($T_{\mathrm{X}}\simeq 14.5\rm{\,keV}$). Model Bl_rc24 has the lowest
temperature ($T_{X}\simeq 11\rm{\,keV}$), which is not surprising given that
this is the model with the lowest collisional energy. Finally for all of the
models there is a remarkable agreement between the measured X-ray
luminosities, as well as with the observational value, which is somewhat
unexpected given the dependency of $L_{X}$ on the squared density.

We now proceed to extract X-ray surface brightness profiles from our fiducial
merger runs. We evaluate these profiles by cutting spatial regions across the
wake behind the X-ray emission peak of the secondary. As in Molnar &
Broadhurst (2015) we identify these regions by visual inspection, with their
spatial location and orientation being changed until the observed spectra best
reproduce previous findings. These regions correspond to the dashed orange
lines shown in the X-ray images of the fiducial models.

Figure 6 shows the measured profiles as a function of the distance $\Delta$
across the wake; these can be compared with those depicted in Figure 8 of
ZYL15 or Figure 4 of Molnar & Broadhurst (2015), after appropriate rescaling.
All of the profiles agree quite well with each other, they share a minimum at
$\Delta\sim-10\rm{\,arcsec}$ and their two peaks are separated by about $\sim
50\rm{\,arcsec}$. The profiles are asymmetric around $\Delta=0$, unlike the
ones presented by ZYL15 and Molnar & Broadhurst (2015). This is easily
explained by a miscentering of the grid used to evaluate the spectra.

Table 4: Main cluster properties for the fiducial merger models.

| B_1 | Bf_rc20 | Bg_rc20 | Bl_rc24 | Observations a
---|---|---|---|---|---
$d_{DM}[\rm{\,kpc}]$ | 704 | 707 | 699 | 707 |  $\sim$ 700
$L_{X}[10^{45}[\rm{\,ergs^{-1}}]$ | 2.86 | 2.75 | 2.84 | 2.15 | 2.19 $\pm$ 0.11
$T_{X}[\rm{\,keV}]$ | 11.8 | 13.2 | 13.5 | 10.9 | 14.5 $\pm$ 1
$d_{X-SZ}[\rm{\,kpc}]$ | 504 | 391 | 414 | 398 |  $\sim$ 600
$\delta T[\mu\kern 2.0pt{}^{\circ}{\kern-2.0ptK}]$ | -915 | -959 | -997 | -662 | -1046 $\pm$ 116
$V_{r}[\rm{\,km\leavevmode\nobreak\ s^{-1}}]$ | 1024 | 956 | 971 | 856 | 586 $\pm$ 96

Notes. Rows from top to bottom: projected distance between the mass centroids,
total X-ray luminosity in the $[0.5-2]\rm{\,keV}$ band, X-ray spectral
temperature of the SE cluster, offset between the SE X-ray peak and the SZ
centroid, central SZ decrement, mean relative radial velocity between the SE
and NW clusters. a: From Table 2 of ZYL15.

The separation between the two maxima and their emission level ($\sim 1.5\cdot
10^{-8}\linebreak\rm{\,counts\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\
cm^{-2}\leavevmode\nobreak\ arcsec^{-2}}$) are in accord with those exhibited
by model B of ZYL15. However for their model the X-ray emission of the minimum
is almost at peak levels, whereas here it is a factor $\sim 2$ lower. This
discrepancy is less pronounced in the spectra of Molnar & Broadhurst (2015),
for which the emission level of the minimum is about $\sim 70\%$ of the
maxima.

The SZ-X-ray offset ($d_{X-SZ}\sim 500\rm{\,kpc}$ ) is substantial in accord
with the corresponding value reported in their Table 2 by ZYL15 for their
fiducial model B, and it is smaller than the value estimated by observations
($\sim 600\rm{\,kpc}$). Following ZYL15 we argue that such a discrepancy is
primarily due to the low angular resolution of the SZ measurements. For the
other merger models the discrepancy becomes larger, as a consequence of the
reduced collision velocity and of the primary’s mass.

The SZ central temperature decrement of fiducial model B is
$\simeq-915\mu\kern 2.0pt{}^{\circ}{\kern-2.0ptK}$, in line with observational
estimates and in accord with the corresponding value reported by ZYL15 for
their model B ($-1130\mu\kern 2.0pt{}^{\circ}{\kern-2.0ptK}$). Models Bf_rc20
and Bg_rc20 have similar decrements, whilst for model Bl_rc24 the decrement is
much smaller ($\simeq-660\mu\kern 2.0pt{}^{\circ}{\kern-2.0ptK}$) and
inconsistent with the measured value. We consider this discrepancy a clear
consequence of the reduced collisional energy of this merging model.

Another observational feature which can be compared against simulation results
is the relative radial velocity $V_{r}$ along the line-of-sight between the NW
and SE clusters. From the redshift distribution of member galaxies, Menanteau
et al. (2012) estimate for the SE cluster a relative radial velocity of
$V_{r}^{s}=598\pm 96\rm{\,km\leavevmode\nobreak\ s^{-1}}$, with respect to the
NW component.

The merging runs presented in this Section, as well as those of ZYL15, do not
incorporate a stellar component in the initial condition setup. Therefore in
these simulations the observed galaxy velocity distribution can only be
contrasted against simulation results if one assumes the velocity distribution
of DM particles to be a proxy of that observed for galaxies. This is supposed
to be a fair assumption as long as the DM is supposed to be collisionless. We
will investigate in Section 3.3 the observational consequences when this
assumption is no longer valid.

For the four fiducial merger models we show in Figure 7 histograms of the
velocity distribution of DM particles, for both NE and SW clusters. The
histograms are normalized to a total of 51 (36) members for the NW (SE)
component, as in Figure 9 of of Menanteau et al. (2012). The Figure also
reports the difference $V^{DM}_{r}$ along the line-of-sight between the mean
radial velocities of the NW and SE clusters. The mean projected velocities of
each component are evaluated by averaging over all of the DM particles which
in the observer plane are within a distance of $400\rm{\,kpc}$ from its mass
centroid.

The velocities are in the range from $V^{DM}_{r}\sim
1,020\rm{\,km\leavevmode\nobreak\ s^{-1}}$ (model B) to $V^{DM}_{r}\sim
850\rm{\,km\leavevmode\nobreak\ s^{-1}}$ (model Bl_rc24) and are significantly
higher than the measured value $V_{r}^{s}$. This tendency has already been
noted by ZYL15. The authors argue that contamination between the different
mass components, which occurs during the collision between the two clusters,
can significantly bias the measurement of radial velocities and cause
underestimation of the relative real cluster velocity.

Figure 6:  X-ray surface brightness profiles across the wake of several
cluster merger models, the extraction region of the profiles being shown in
the previous images as the orange dashed line. Distance is measured in arc sec
along the extraction path, model B_1 is the corresponding B model of ZYL15 and
is used here as a reference model.

The two velocity histograms of fiducial model B can be directly compared with
the corresponding ones shown in Figure 9 (bottom panel) of ZYL15, the
construction procedure being the same. From the two Figures it appears that
the velocity distributions extracted from the two simulations are in good
agreement. In particular ZYL15 report $V^{DM}_{r}\sim
910\rm{\,km\leavevmode\nobreak\ s^{-1}}$ for their model B.

To summarize, we use the results of this Section to argue that the observed
twin-tailed X-ray structure seen in the merging cluster El Gordo can be
reproduced in off-axis merging simulations only for a narrow interval of
values of the initial collision parameters. Specifically, the observational
morphological constraints are satisfied for $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}$. The mass of the
primary is allowed to vary in the range $10^{15}\rm{\,M_{\odot}}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}M^{(1)}_{200}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2.5\cdot 10^{15}\rm{\,M_{\odot}}$,
with the mass ratio $q$ adjusted to satisfy $M^{(2)}_{200}\simeq 6.5\cdot
10^{14}\rm{\,M_{\odot}}$. The lower part of the mass interval is favored by
the latest mass estimates.

The fiducial merger models which have been found here to agree with the
observed X-ray structures also satisfy a number of observational constraints.
In particular models Bf_rc20 and Bg_rc20, which have a primary mass of
$M^{(1)}_{200}=1.6\cdot 10^{15}\rm{\,M_{\odot}}$, show an overall better
agreement with data.

Figure 7: Histograms of the DM line-of-sight velocity distribution for the NW
and SE clusters of the same merger models shown in Figure 6. The NW (SE)
histograms are normalized to a total number of 51 (36) members, as in Figure 9
of ZYL15. For each model in the Figure is reported the relative mean radial
velocity between the SE and NW cluster components.

###  3.2 A returning scenario ?

One of the most challenging features to be reproduced in the simulations of
the El Gordo cluster is the observed offset between the X-ray emission and DM
density peak of the SE component, with the latter being closer to the merging
center of mass (see, for example, Figure 1 of Ng et al. 2015). This
observational aspect is clearly at variance with what is expected in a
outgoing, post-pericenter, scenario in which, after the collision, the X-ray
emission peak is expected to trail the DM mass centroid.

To solve this inconsistency Ng et al. (2015) suggested that the El Gordo
cluster is presently observed at a later stage of the merger, in a returning
phase with the DM secondary having reached the apocenter and now falling back
onto the primary. The authors conclude that a returning scenario is favored by
comparing available data against a large set of Monte-Carlo simulations with
different input parameters. They construct different merger models by
following the evolution of the two DM subclusters in a head-on scenario.

In this Section we perform a suite of hydrodynamical simulations to
investigate the viability of a returning scenario as a possible solution to
the offset problem. To consistently compare our simulations with previous
findings (Ng et al. 2015) we will consider head-on collisions and set
accordingly $P=0$, the procedures for setting up the initial conditions being
that given in Section 2.

We construct our set of simulations by considering a range of collision
parameters. Specifically, for the mass of the primary we choose
$M^{(1)}_{200}=1.3\times 10^{15}\rm{\,M_{\odot}}$, as in Table 1 of Ng et al.
(2015). The mass of the secondary is set by those authors to
$M^{(2)}_{200}=7.6\times 10^{14}\rm{\,M_{\odot}}$, or $q\simeq 1.7$. Here we
instead consider the following values for the merging mass ratio:
$q=2,\leavevmode\nobreak\ 4,\leavevmode\nobreak\ 10$.

For a given value of $q$ we consider three different values of the collision
velocity: $V=1,000,\leavevmode\nobreak\ 2,000$ and
$3,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$, so that our set of simulations
consists of nine head-on merger collisions. The IDs of the runs are given in
Table 5, these are defined according to the mass ratio $q$ and the initial
collision velocity $V$. We explicitly decided to investigate a wide range of
merging parameter space, in order to assess for which initial conditions head-
on encounters are best able to reproduce the observed twin-tailed X-ray
morphology.

As the collision parameter $P\rightarrow 0$ the asymmetries in the X-ray
emissivity maps of the merging clusters tend to disappear and the collision
becomes axially symmetric. Previous head-on merging simulations of the El
Gordo cluster have been conducted by Donnert (2014) and ZYL15 (their model A).
As can be seen from Table 1, these simulations were not perfectly head-on but
were performed with an initially small impact parameter ($P=300\rm{\,kpc}$).

The simulations of Donnert (2014) were able to satisfy a number of
observational constraints but failed to reproduce the twin-tailed X-ray
morphology, the X-ray peak being followed by only one tail elongated along the
collision axis. Model A of ZYL15 is a very energetic
($V=3,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$), almost head-on collision. As
discussed in the previous Section, their fiducial model A corresponds to the
bottom mid panel of Figure 1 and the twin-tailed morphology is absent in the
highly symmetric X-ray structure. A small asymmetric wake-like feature is
instead present in the bottom right panel of Figure 1 (Figure 1c of ZYL15),
which the authors interpret as a consequence of a non-zero impact parameter.

In the latter case it is worth noting the different values found in the two
simulations for the projected separation between the mass centroids. ZYL15
report $d_{DM}\simeq 600\rm{\,kpc}$, while we found here $d_{DM}\simeq
720\rm{\,kpc}$. We attribute this discrepancy to the very large value of the
relative velocity ($V_{rel}\simeq 4,700\rm{\,km\leavevmode\nobreak\ s^{-1}}$)
between the two clusters at the evolution time $t=0.09$ Gyr, when the maps are
extracted, and the time resolution of our simulations, given by the timestep
$\Delta t_{0}=1/100$ Gyr. The current epoch is identified at run time when
$d_{DM}\geq 700\rm{\,kpc}$, and the snapshot of Figure 1 (bottom right)
corresponds to the simulation step $n=164$. The previous output is saved at
$n=162$ , when $d_{DM}\simeq 630\rm{\,kpc}$.

Table 5:  IDs and initial parameters of the head-on merging simulations of
Section 3.2 (a)

$V[\rm{\,km\leavevmode\nobreak\ s^{-1}}]$ | $q=M_{1}/M_{2}$ |
---|---|---
| 2 | 4 | 10 |
$1,000$ | R01V10 | R02V10 | R03V10 |
$2,000$ | R01V20 | R02V20 | R03V20 |
$3,000$ | R01V30 | R02V30 | R03V30 |

Notes. a: The meaning of the parameters $q$ and $V$ is the same as in Table 1.
For all the simulations we set the mass of the primary to
$M^{(1)}_{200}=1.3\cdot 10^{15}\rm{\,M_{\odot}}$ and the initial gas mass
fraction to $f_{g}=0.1$ for both of the clusters.

To summarize, no twin-tailed X-ray structure has been reproduced so far in
head-on merging simulations of the El Gordo cluster. We thus now discuss
results extracted from our set of head-on merging simulations. For each of the
merging runs of Table 5 we show in Figure 8 the corresponding X-ray emission
maps in the outgoing scenario. The collisions are seen face-on and the
secondary moves from right to left. The time $t$ displayed in each panel is
when for the first time $d_{DM}\simeq 700\rm{\,kpc}$ after the pericenter
passage; here $t=0$ is when the simulations start.

Figure 8:  Surface brightness maps for the head-on mergers described in
Section 3.2. As in the previous maps time is in Gyr, with the pericenter
passage defined at $t=0$ and the box size is $1.6\rm{\,Mpc}$. In each panel is
indicated the collision model whose initial merging parameters are given in
Table 5. Each row of panels refers to simulations with the same value of the
mass ratio $q$, while each column is for simulations having the same value of
$V$. This is reported in $\rm{\,km\leavevmode\nobreak\ s^{-1}}$. The maps are
extracted at output times when the distance $d_{DM}$ between the mass
centroids of the two clusters is approximately $d_{DM}\simeq 700\rm{\,kpc}$.
As in Figure 1, the crosses indicate the projected spatial locations of the
mass (green) and X-ray surface brightness (red) centroids.

From Figure 8 it can be seen that all of the simulations shown are
axisymmetric around the collision axis. Moreover, not unexpectedly, the X-ray
emission peaks are always trailing the DM centroids. The only exception to
this is the run R02V02, but its initial conditions are a bit peculiar (see
later). Finally, in all of the cases, the core of the primary does not survive
the encounter with the secondary, with the ICM behind the secondary being
shaped by its motion through the main cluster. For the R01 runs ($q=2$), the
ICM morphology is arc-shaped around the secondary’s core, in accord with
previous findings (Ricker & Sarazin 2001), while for higher mass ratios the
post-pericenter shape of the ICM depends on the initial velocity $V$.

Figure 9:  As in Figure 8, but here the output times are chosen so that
$d_{DM}\simeq 700\rm{\,kpc}$ after the apocenter passage.

The only merging simulations which show X-ray structures whose shapes can be
considered to exhibit some degree of similarity with the requested twin-tailed
morphology are the R01 runs ($q=2$). In fact, the initial conditions of the
R01V30 run are very close to those of model A in Table 1, while those of
Donnert (2014) can be considered intermediate between R01V30 and R01V20. In
the other mergers with $q>2$, the smaller is the mass of the secondary, the
less is the impact on the ICM of the primary. A significant X-ray structure
behind the secondary is still visible if the decrease in mass of the secondary
is accompanied by a corresponding reduction in the collision velocity (runs
R02V10 and R03V10).

To assess the viability of the returning scenario, Figure 9 shows X-ray maps
extracted from the same runs as in Figure 8, but now in a post-apocenter
phase. The criterion to identify the observer epoch is still that of having
$d_{DM}\simeq 700\rm{\,kpc}$, but this when the secondary has reached the
apocenter and is now returning towards the primary.

The time difference $\Delta t_{orb}$ between the epochs shown in Figure 9 and
those in the corresponding panels of Figure 8 gives, for a specific model, the
time spent by the secondary between the outgoing phase and the returning one.
The time $\Delta t_{orb}$ can assume a wide range of values: from $\Delta
t_{orb}\sim 0.02$ (R02V02) up to $\Delta t_{orb}\sim 4$ Gyr (R03V30).

The time $\Delta t_{orb}$ is a function of both $q$ and $V$; for a given $V$,
it can be seen that it increases as $q$ gets higher. Such behavior is clearly
a consequence of the reduced exchange of energy between the two clusters
during the collision. The merging run R02V02 has a very small value of $\Delta
t_{orb}$ because its initial velocity ($V=2,000\rm{\,km\leavevmode\nobreak\
s^{-1}}$) is such that when the secondary reaches the orbit’s apocenter it is
also when $d_{DM}\simeq 700\rm{\,kpc}$; this is the reason why the X-ray
emission and DM centroids of the secondary are so close in both of these two
phases.

To compare our findings with the proposed model of Ng et al. (2015) we must
choose from Figure 9 those merging runs for which the value of the mass ratio
is closer to the value ($q\sim 1.7$) adopted by them. These conditions are
satisfied by the runs R01V10, R01V20 and R01V30, for which $\Delta t_{orb}\sim
0.8,0.4$ and $\sim 1$ Gyr, respectively. The latter model is however ruled out
by the timescale $\Delta t_{orb}=0.91-0.46\sim 0.5$ Gyr favored by Ng et al.
(2015). Their Monte-Carlo simulations include neither dynamical friction nor
tidal forces, so that a comparison between their favored value of $\Delta
t_{orb}$ with the values we found here is not entirely consistent.
Nonetheless, in a previous paper (see Table 2 of Valdarnini & Sarazin 2021) we
demonstrated that the impact of dynamical friction and tidal forces becomes
significant for off-axis mergers when the mass ratio gets higher ($q\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}5$). We then use these findings to
argue that the constraint on $\Delta t_{orb}$ of Ng et al. (2015) can be used
to rule out model R01V30.

Table 6: Initial merger parameters and IDs of the merger simulations of Figure
10 and 11. a

Model: $\\{M_{1},r_{s},\leavevmode\nobreak\ q,\leavevmode\nobreak\ P,\leavevmode\nobreak\ V,\leavevmode\nobreak\ \zeta=r_{s}/r_{c}\\}$ | BCGs | $\sigma_{DM}/m_{X}[\rm{\,cm^{2}}{gr^{-1}}]$ | IDs of the merger simulations
---|---|---|---
Bf$\\_$rc20: $\\{1.6\cdot 10^{15},0.696,\leavevmode\nobreak\ 2.32,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 2,500,\leavevmode\nobreak\ 3.44\\}$ | $\surd$ | 0 | DBf$\\_$rc20
Bl$\\_$rc24: $\\{1.\cdot 10^{15},0.574,\leavevmode\nobreak\ 1.54,\leavevmode\nobreak\ 600,\leavevmode\nobreak\ 2,000,\leavevmode\nobreak\ 2.39\\}$ | $\surd$ | 0 | DBl$\\_$rc24
Bf$\\_$rc20 | $\surd$ | 1/2/5 | XDBf$\\_$rc20
Bl$\\_$rc24 | $\surd$ | 1/2/5 | XDBl$\\_$rc24

Notes. a Columns from left to right: collision parameters of the mergers,
cluster halos with initially a BCG mass component, value of the self-
interacting DM cross section per unit mass, ID of the merger simulations. For
the simulations with ${\sigma_{DM}}/{m_{X}}>0$, the label X generically refers
to merger runs with different values of ${\sigma_{DM}}/{m_{X}}$.

This leaves us with models R01V20 and R01V10. Note that for model R01V20 the
DM centroid of the secondary is now closer to the center of mass of the system
than its X-ray emission peak, the distance between the two centroids being
approximately $\sim 100\rm{\,kpc}$. However neither of these two merger models
is able to reproduce the X-ray structure that we see in the El Gordo clusters.

In fact these models, as well as all of the others, are morphologically
inconsistent with the observed twin-tailed X-ray morphology. We argue that
this behavior can be interpreted as follows. The primary’s ICM undergoes a
strong increase in temperature and luminosity as the two cluster cores
approach the pericenter, when their relative separation has a minimum.
Subsequently the two clusters move away from each other and the shock-heated
gas, which is no longer in hydrostatic equilibrium, cools adiabatically and
its temperature and density drop significantly. Because of the X-ray
emissivity dependence on the squared density, this in turn implies a
substantial reduction in the X-ray brightness of the gas structures created
during the collision.

This behavior is in line with previous findings (Ricker & Sarazin 2001; Poole
et al. 2006) and in particular with ZYL15, who argue that the two wings seen
at $t=0.09$ after the pericenter passage in their model of Figure 1c are
short-lived. However the X-ray brightness maps presented by Donnert (2014,
Figure 5 ) in his simulation study of El Gordo present some degree of
persistency even at late times ($t\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}2$ Gyr). The origin of this
discrepancy it is not clear; we suggest that since the Donnert (2014)
simulations were made using standard SPH, shortcomings of the code in the
modeling of hydrodynamical instabilities might have played a role.

The results presented here are valid for head-on collisions, nonetheless it
should be emphasized that off-axis mergers as well are clearly inconsistent
with a returning scenario. This can be clearly seen by realizing that for
these merger models the time difference $\Delta t_{orb}$ between the pre- and
post-apocenter epochs, at which must be satisfied the constraint $d_{DM}\simeq
700\rm{\,kpc}$, are significantly higher than for mergers with orbits having
zero angular momentum. For the fiducial models of Section 3.1 it is found
$\Delta t_{orb}\sim 1-2$ Gyr, whereas we have seen that the post-collision
X-ray structures are fast evolving and their lifetime lies in the range $\sim
0.1-0.2$ Gyr.

To summarize, the results of this Section strongly support the outgoing
scenario to explain the observed features of the merging cluster El Gordo. In
the next Section we will investigate alternative possibilities to reproduce
the observed offsets between the different mass components of the El Gordo
cluster.

###  3.3 Mergers with BCGs and self-interacting DM

One of the most interesting aspects of the merging cluster El Gordo is the
measured offset between the projected positions of the X-ray emission and DM
centroid of the SE cluster ($\sim 100\rm{\,kpc}$). As already said in the
Introduction, the relative position of the X-ray peak with respect to the mass
distribution is in contrast with what is expected from dissipative arguments.
Moreover, non-zero offsets are also measured (Kim et al. 2021) between the BCG
to DM centroid ($\sim 60\rm{\,kpc}$) and the BCG to X-ray peak ($\sim
70\rm{\,kpc}$) centroid.

According to the ‘central galaxy paradigm’ (van den Bosch et al. 2005) BCGs
are expected to reside at the center of the DM halo, thus tracing the bottom
of the potential. In their simulation study, Martel et al. (2014) concluded
that the paradigm is statistically satisfied, but that exceptions are common.
From a large survey ($\sim 10,000$) Zitrin et al. (2012) found a relatively
small mean offset ($\sim 15\rm{\,kpc}$) between the BCG and DM cluster
centroids. More recently Depropis et al. (2021) showed that a significant
fraction of BCGs are offset from the X-ray centroid, but are still aligned
with the mass distribution.

In order to investigate the impact of a massive merging cluster like El Gordo
on the ‘central galaxy paradigm’, i.e. that the BCGs are located at the
cluster center, we resimulate two of the fiducial models of Section 3.1. We
choose fiducial models Bf_rc20 and Bl_rc24, whose simulation results produced
overall the best fitting to the data and which have a mass of the primary that
brackets the observational range.

Our new models will be labeled as DBf_rc20 and DBl_rc24, the initial condition
setup being the same as of their respective models in Section 3.1, but each
cluster now has a stellar mass component added to its initial mass
distribution. A numerical realization of the star density profiles which are
supposed to represent the BCGs is described in Section 2.2.3. The stellar
masses of the primary and secondary cluster that we obtain from Equation (14)
are $M^{(1)}_{\star}=2.2\cdot 10^{12}\rm{\,M_{\odot}}$ and
$M^{(2)}_{\star}=1.6\cdot 10^{12}\rm{\,M_{\odot}}$, respectively. Table 6
lists the main merger parameters of the two models.

Figure 10 shows the mock X-ray images extracted from the two runs at the
simulation time identified as the present epoch. It is easily seen that adding
a BCG stellar component to the two clusters does not modify in any way the
X-ray morphology of the mergers. The two X-ray maps are very similar to the
corresponding maps shown in the Bf_rc20 and Bl_rc24 panels of Figure 2 and
Figure 4, respectively.

Figure 10:  X-ray surface brightness images extracted from simulations of
models DBf_rc20 and DBl_rc24. The merging parameters are the same of models
Bf_rc20 and Bl_rc24, respectively, but the runs were performed with initial
conditions including for the two clusters a stellar mass component describing
a brightest cluster galaxy (see Section 2.2.3). In each panel the open orange
stars mark the projected spatial location of the mass centroids of the star
particles representing the BCGs. The value of $V^{s}_{r}$ is relative mean
radial velocity along the line-of-sight between the SE and NW BCG components.
The meaning of the other symbols and the viewing direction is the same of the
maps shown in Figure 1.

About the dynamics of the BCGs two important results emerge from Figure 10.
The first concerns the measured velocity offsets $V_{r}^{s}$ along the line-
of-sight between the two BCGs of each merger model. As can be seen from Figure
7, these values are relatively speaking very close to the velocities $V_{r}$
extracted from the DM particle velocities of the corresponding model without a
star component. This confirms that in a collisionless DM merging scenario the
BCG star velocities can be assumed to be a fair proxy for the DM particle
velocities.

The other important result is that the mass centroids of the BCGs appear to be
strictly aligned with their corresponding DM mass centroids. This indicates
that the off-center positions of the BCGs that we observe today were not
generated by the merging process during the cluster collision. According to
Martel et al. (2014) the spatial offset between the BCG mass centroid and that
of its parent cluster is produced in major cluster mergers, as a result of the
collisions. If at $z=0$ the cluster has reached an equilibrium configuration,
the BCG will have settled at the bottom of the potential and its off-center
location will be much reduced or even absent. We will return in the
conclusions to the likelihood of such a scenario for the merging cluster El
Gordo.

An alternative scenario which could explain the observed spatial offsets
between the different mass components is that of having DM which is self-
interacting, these signatures being a characteristic of SIDM in major mergers.
These expected properties have led many authors to investigate cluster
collisions in an SIDM framework (Markevitch et al. 2004; Kahlhoefer et al.
2014; Kim et al. 2017; Robertson et al. 2017, 2019; ZuHone et al. 2019; Fisher
et al. 2022).

We investigate here the possible outcomes of SIDM cluster mergers, using the
implementation of the SIDM modeling described in Section 2.3. We performed
SIDM merging simulations for each of the two merger models of Figure 10, and
we generically refer to these models as XDBf_rc20 and XDBl_rc24, respectively.
Three SIDM simulations were carried out for each model, with the self-
interacting DM cross section assuming the following values:
${\sigma_{DM}}/{m_{X}}=1,\leavevmode\nobreak\ 2$ and
$5\rm{\,cm^{2}}{gr^{-1}}$, respectively. In Table 6, we list all models with
their respective initial merger parameters.

The mock X-ray maps extracted from the SIDM runs are presented in Figure 11,
with the top panels being for the model XDBf_rc20 and bottom panels for model
XDBl_rc24. We also considered SIDM simulations with
${\sigma_{DM}}/{m_{X}}=0.1\rm{\,cm^{2}}{gr^{-1}}$, but the results are almost
indistinguishable from the standard CDM mergers and are not shown here.

A first important conclusion to be drawn from the maps of Figure 11 concerns
the behavior of the X-ray gas morphology. In an SIDM scenario the interaction
of DM leads to an exchange of energy between the two clusters during the
collision, which implies the development of shallower DM potential wells for
the two DM halos (Kaplinghat et al. 2016; Kim et al. 2017; ZuHone et al.
2019). This can be easily seen by contrasting the contour levels of the
projected mass density: those displayed in the map of Figure 11 from the
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ run of model XDBf_rc20 appear
much rounder than those exhibited in Figure 2 by their counterparts from the
standard CDM merger Bf_rc20.

The presence of shallower DM potential wells leads in turn to a reduced
resiliency of the post-pericenter gas structures, which can now more easily
escape from the potential wells of the original halos. As can be seen from
Figure 11, this effect is present in all of the runs, but increases as
${\sigma_{DM}}/{m_{X}}$ gets higher and when the mass of the primary gets
lower. In fact, the bottom right panel of Figure 11 shows that for
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ model XDBl_rc24 has entirely
lost its original X-ray morphology. A complementary role in the process
causing the disappearance of the post-pericenter gas structures is played by
the elapsed time $t$ between the pericenter passage and the observer epoch.
This time is much larger ($\sim 0.2$ Gyr) than in the collisionless DM runs,
and thus the X-ray structures generated during the collision are much weaker
when observed or even absent (see Section 3.2).

Figure 11:  X-ray surface brightness maps extracted at the present epoch $t$
from two sets of SIDM merging runs. The simulations of the top (bottom) panels
have the same initial conditions of model DBf_rc20 (DBl_rc24) of Figure 10.
From left to right each row of panels corresponds to SIDM simulations
performed for three different values of the DM scattering cross-section:
${\sigma_{DM}}/{m_{X}}=1,\leavevmode\nobreak\ 2$ and
$5\rm{\,cm^{2}}{gr^{-1}}$, respectively. In each panel the distance
$d_{\rm{X-DM}}$ indicates the value in kpc of the projected distance between
the X-ray emission peak and the DM mass centroid, $d_{\rm{BCG-X}}$ that
between the mass centroid of the BCG galaxy and the X-ray emission peak, and
finally $d_{\rm{BCG-DM}}$ is the distance between the BCG and DM mass
centroids. All of the centroids refer to the SE cluster. The meaning of the
other symbols is the same as in Figure 10. The thresholds and the spacing of
the contour levels are the same as those in Figure 1. In each panel the value
of $n$ is the simulation step at which the map is extracted and $t=0$ at the
pericenter passage. The observer epoch is defined as in the other Figures.

This effect can be clearly seen in the panels of Figure 11, where each panel
reports the time $t$ since passing pericenter. For the XDBf_rc20 simulations
(top panels) the DM scattering cross-section increases from
${\sigma_{DM}}/{m_{X}}=1\rm{\,cm^{2}}{gr^{-1}}$ up to
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$, while the time $t$ has a
corresponding increase from $t=0.15$ Gyr up to $t=0.24$ Gyr. For XDBl_rc24 the
effect is even more pronounced, with $t$ now ranging from $t=0.19$ Gyr up to
$t=0.45$ Gyr. To better illustrate the difference from the standard CDM runs
of Figure 10, we also report in each panel the simulation step $n$. From this
one can recover the simulation time $t_{n}=n\Delta t_{0}$, at which the maps
were extracted.

We attribute this increasing dependency of the elapsed time $t$ on the DM
scattering cross-section ${\sigma_{DM}}/{m_{X}}$ to the mutual transfer of
energy which occurs between the two DM halos during the collision. In the off-
axis mergers DBf_rc20 and DBl_rc24 the DM halos remain almost unperturbed, but
in the SIDM runs the exchange of energy can be significant when
${\sigma_{DM}}/{m_{X}}\sim 5\rm{\,cm^{2}}{gr^{-1}}$. As a consequence of these
dissipative processes, the post-pericenter cluster bulk velocities are now
strongly reduced with respect the standard CDM model. As will be shown below,
this in turn will also impact on the BCG bulk velocities.

Another important feature which should be present when observing merging
clusters in an SIDM scenario are the expected offsets between the different
mass components. Because of the dissipative behavior of the DM, during the
merger the collisionless galaxy component is expected to decouple from its
parent DM halo (Kim et al. 2017). Moreover, the position of the X-ray peak in
relation to that of the DM centroid is also expected to differ in a
significant way from that found in the corresponding standard CDM merger
(Tulin & Yu 2018).

In the mock X-ray maps of Figure 11, for each merging run we indicate in the
corresponding panel the projected spatial positions of the different
centroids. As in the X-ray maps of Figure 10, the red crosses indicate the
positions of the X-ray peaks, while the green crosses are for the DM
centroids. The orange stars refer to the centroids of the BCGs.

From the merging simulations of Figure 11, there are several noteworthy
features which are clearly recognizable by looking at the relative positions
of the different centroids. As expected, the positions of the BCG centroids
are now decoupled from those of their parent DM halos. Moreover, for the SE
cluster the position of the X-ray peak is highly significant. For the standard
CDM mergers of Figure 10 the DM peak is located farther away from the
barycentre than the X-ray centroid, and the relative offset is negative. For
the SIDM merging simulations of Figure 11 the relative separation between the
two peaks initially decreases as ${\sigma_{DM}}/{m_{X}}$ increases, and
becomes positive for ${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$. For
these runs the X-ray centroid is no longer lagging behind the mass centroid
but is now clearly ahead of it.

Table 7: IDs and initial merger parameters of the two SIDM merging simulations
of Figure 12 .a

Model | $M^{(1)}_{\star}[\rm{\,M_{\odot}}]$ | $M^{(2)}_{\star}[\rm{\,M_{\odot}}]$ | $N^{(1)}_{\star}$ | $\varepsilon_{\star}[\rm{\,kpc}]$ | $r_{c}^{\rm{prm}}[\rm{\,kpc}]$ | $\zeta$ | $\sigma_{DM}/m_{X}[\rm{\,cm^{2}}{gr^{-1}}]$ | $(f_{g1},f_{g2})$
---|---|---|---|---|---|---|---|---
XDBf_sa | $2.2\times 10^{12}$ | $1.6\times 10^{12}$ | $16,785$ | 9.5 | 406 | 1.7 | 5 | (0.12,0.12)
XDBf_sb | $2.2\times 10^{12}$ | $1.6\times 10^{12}$ | $16,785$ | 9.5 | 290 | 2.4 | 4 | (0.12,0.14)

Notes. a Columns from left to right: ID of the merging model, stellar mass of
the BCG of the primary, the same mass but for the secondary, number of star
particles for the primary, gravitational softening length of the star
particles, gas core radius of the primary, dimensionless parameter
$\zeta=r_{s}/r_{c}$, self-interacting DM cross section per unit mass, primary
and secondary cluster gas mass fractions $f_{g}$ at $r_{200}$. For the two
SIDM merger models the collision parameters are those of model Bf in Table 1:
$\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
V,\leavevmode\nobreak\ P\\}=\\{1.6\cdot
10^{15}\rm{\,M_{\odot}},\leavevmode\nobreak\ 2.32,\leavevmode\nobreak\
2,500\rm{\,km\leavevmode\nobreak\ s^{-1}},\leavevmode\nobreak\
600\rm{\,kpc}\\}$ .

For the same value of ${\sigma_{DM}}/{m_{X}}$, the relative separations
between the different centroids should be smaller in model XDBl_rc24 than in
model XDBf_rc20, the collision velocity being smaller
($V=2,000\rm{\,km\leavevmode\nobreak\ s^{-1}}$) in the former than in the
latter ($V=2,500\rm{\,km\leavevmode\nobreak\ s^{-1}}$). Leaving aside the
${\sigma_{DM}}/{m_{X}}=1\rm{\,cm^{2}}{gr^{-1}}$ runs, for which this
dependency is too weak to be detected, the two merging simulations with
${\sigma_{DM}}/{m_{X}}=2\rm{\,cm^{2}}{gr^{-1}}$ (middle panels of Figure 11)
effectively reproduce this behavior.

However for the merging runs with
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ the separations between the
centroids of model XDBl_rc24 are much larger that those of model XDBf_rc20. We
argue that this discrepancy is due to an observational effect, with the
elapsed time $t$ being much larger in model XDBl_rc24 ($t\sim 0.45$ Gyr) than
in model XDBf_rc20 ($t\sim 0.24$ Gyr).

The same time difference is responsible of another noteworthy feature of model
XDBl_rc24 with ${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$. From the
bottom right panel of Figure 11 it can be seen that this run exhibits negative
BCG-DM offsets, i.e. the DM peaks are more distant from the barycentre than
the BCGs.

As the clusters collide, the DM particles will lose momentum through
scatterings. The resultant drag force will lead to a deceleration of the DM
halos and in turn to positive offsets between the DM and galaxy components
(Markevitch et al. 2004; Kim et al. 2017; Robertson et al. 2017; Fisher et al.
2021, 2022). However this description applies only in proximity of the first
pericenter passage. At later times the gravitational pull of a DM halo will
act on the corresponding galaxy component as a restoring force, which will
thus decelerate and reduce its bulk motion (Harvey et al. 2014).

According to the relative strength of the DM cross-section and the elapsed
time t, from initially positive the BCG-DM offsets can then become negative,
as seen in the bottom right panel of Figure 11 for model XDBl_rc24. These
behaviors are in line with what is expected in a merging SIDM scenario, being
in general non-linear the dependency on ${\sigma_{DM}}/{m_{X}}$ of the
relative distance between the DM centroid and that of the other mass
components.

A significant conclusion that can be drawn from these findings is that model
XDBf_rc20 with ${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ is the only
merging simulation which is able produce relative separations within the
observational range ($\sim 50-100\rm{\,kpc}$). The validity of such a
conclusion being clearly dependent on the statistical significance of the size
of the observed offsets for the El Gordo cluster, we will come back later to
this issue.

As a consequence of this conclusion the SIDM model adopted here to simulate
the merging cluster El Gordo has no free parameters, the only admitted values
for the DM scattering cross-section being in the range
${\sigma_{DM}}/{m_{X}}\sim 4-5\rm{\,cm^{2}}{gr^{-1}}$.

The reason for this limitation in SIDM merging simulations of the cluster El
Gordo is that the relative separations between the centroids of the different
mass components are locked to the distance $d_{DM}$ between the DM mass
centroids, which fixes the present epoch through the constraint $d_{DM}\sim
700\rm{\,kpc}$.

This is made clear by the middle panels of Figure 11: where for
${\sigma_{DM}}/{m_{X}}=2\rm{\,cm^{2}}{gr^{-1}}$ both models exhibit a relative
separation between the different centroids of about $\sim 50\rm{\,kpc}$, at
approximately the same epoch $t$. However for
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ only model XDBf_rc20 is able
to reproduce the requested observational range for the relative separations,
while model XDBl_rc24 completely fails to match the data. This behavior as
${\sigma_{DM}}/{m_{X}}\rightarrow 5\rm{\,cm^{2}}{gr^{-1}}$ is due to the
relative impact of DM collisions, which is much higher in model XDBl_rc24 than
in XDBf_rc20. As previously explained, for model XDBl_rc24 this value in turn
implies lower post-pericenter DM bulk velocities and a much higher elapsed
time $t$. Note that these conclusions also work in the other direction: for
SIDM merging simulations of the El Gordo cluster, the best agreement with data
is obtained only for a primary mass of about $M^{(1)}_{200}\sim 1.6\cdot
10^{15}\rm{\,M_{\odot}}$.

Motivated by these results, we now turn our attention to the SIDM merging
simulation of model XDBf_rc20 with
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$, which from the top right
panel of Figure 11 appears as the only viable candidate able to reproduce many
observational properties of the cluster El Gordo. From the contour levels of
the surface mass density (see the caption of Figure 1), we estimate for the SE
cluster of this model a peak value of the projected DM mass surface density of
about $\Sigma_{DM}\sim 0.3\rm{\,gr{\,cm^{-2}}}$. The cluster’s scattering
depth is then $\tau=\Sigma_{DM}{\sigma_{DM}}/{m_{X}}\sim 1.5$, which is
marginally above the optically thick limit (Markevitch et al. 2004). However
this is a peak value. By using an average, one would obtain a lower value for
$\tau$.

As already outlined, for this merging run the DM potential wells of the two
halos are strongly perturbed by the presence of SIDM. The top right panels of
Figure 11 demonstrate that with respect the standard CDM merging simulation
DBf_rc20, the impact of SIDM on the post-pericenter X-ray structures is
significant when ${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$.

We are then forced to reconsider the original initial condition parameters of
the two baryonic halos of model Bf_rc20, in order to reproduce in a
${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$ SIDM merging simulation the
observed twin-tailed X-ray morphology of the El Gordo cluster. The initial
collision parameters of this simulation are those of model Bf in Table 1:
$\\{M^{(1)}_{200},\leavevmode\nobreak\ q,\leavevmode\nobreak\
V,\leavevmode\nobreak\ P\\}=\\{1.6\cdot
10^{15}\rm{\,M_{\odot}},\leavevmode\nobreak\ 2.32,\leavevmode\nobreak\
2,500\rm{\,km\leavevmode\nobreak\ s^{-1}},\leavevmode\nobreak\
600\rm{\,kpc}\\}$.

Figure 12: X-ray images extracted at the observer epoch from two SIDM merging
simulations, the initial conditions of the two runs are given in Table 7. The
meaning of the distances shown in each panel is the same of the ones reported
in Figure 11, with the exception of $d^{NW}_{\rm{SZ-DM}}$ which is the
projected distance between the SZ peak and the DM mass centroid of the NW
cluster. As in Figure 10, $V^{s}_{r}$ is line-of-sight relative mean radial
velocity between the two BCGs. The thresholds and the spacing of the contour
levels are the same as in Figure 1. The filled circles indicate the peak
locations from several observations, as taken from Figure 6 of Kim et al.
(2021). Their spatial positions have been normalized to the relative distance
from the mass centroids (see text). The color coding of the circles is the
same of the associated crosses, which indicate the projected positions of the
corresponding centroids as extracted from the simulations.

For this set of initial collision parameters we performed SIDM merger
simulations for runs with both ${\sigma_{DM}}/{m_{X}}=5\rm{\,cm^{2}}{gr^{-1}}$
and ${\sigma_{DM}}/{m_{X}}=4\rm{\,cm^{2}}{gr^{-1}}$. We found that in order to
produce realistic X-ray morphologies with these simulations, initial
conditions with larger cluster gas fractions and gas core radii of the primary