© The Authors 2024

1. Introduction

Galaxy clusters are the largest virialized structures known in the Universe, and according to the hierarchical scenario, their formation proceeds from the bottom up, through the sequential accretion and merging of smaller structures. They have measured virial masses in the range M ∼ 10¹⁴ − 10¹⁵ M_⊙. Observations show that about ∼80% of their mass content consists of dark matter (DM), with the remaining ∼20% in the form of baryons. The bulk of these baryons (∼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 dissipate the infall kinetic energy and the gas is shock-heated to high temperatures. Because of the masses involved, at equilibrium these temperatures will be as high as T ∼ 10⁷ − 10⁸ 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 be used to study the physics of the ICM and in turn derive cluster properties with which to constrain the currently accepted Lambda cold dark matter (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 for probing 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. Mergers are considered major when the mass ratio of the primary to the secondary is close to unity, and minor whenever the mass of the secondary is ≲30 − 40% of the primary mass. X-ray surface brightness maps of minor mergers 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 are among the most energetic events since the Big Bang, with collision energies in the range ∼10⁶³ − 10⁶⁴ ergs⁻¹ (Sarazin 2002). The energy involved in these collisions renders these objects unique laboratories for probing the collisional properties of their components, as well as the likelihood of such events within the ΛCDM cosmological model (Lee & Komatsu 2010; Lage & Farrar 2015; Asencio et al. 2021).

An important effect that will occur during a major merger is the expected offset between the dissipative ICM and the collisionless components: galaxies and DM. Additionally, self-interacting dark matter (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), the first observed to show 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 nonlinearity 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 that 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 the effects of mergers in general (Roettiger et al. 1996; Ricker & Sarazin 2001; Ritchie & Thomas 2002; Poole et al. 2006, 2008; McCarthy et al. 2007; Mitchell et al. 2009; ZuHone et al. 2009) or 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 aimed at studying a particular merging system to be performed and therefore contrast the simulation results with a given observation. Some examples of such specifically designed simulations include those of the Bullet Cluster (Springel & Farrar 2007; Mastropietro & Burkert 2008; Lage & Farrar 2014), the cluster ACT-CL J0102-4915 (Donnert 2014; Molnar & Broadhurst 2015; Zhang et al. 2015, 2018), 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, 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 on the order of ∼2 × 10¹⁵ M_⊙, with a galaxy velocity dispersion σ_gal ∼ 1300 km s⁻¹ and an integrated temperature in the range T_X ≃ 15 keV. Similar mass estimates have also been obtained from weak lensing (WL; 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 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 ΛCDM model (Menanteau et al. 2012; Jee et al. 2014).

The cluster exhibits a binary mass structure (Jee et al. 2014), with its two subclusters denominated the northwestern (NW) and southeastern (SE), respectively (see, for example, Fig. 1 of Ng et al. 2015) because of their positions. WL analysis (Jee et al. 2014) shows a projected separation between the mass centroids of the two components of approximately d ∼ 700 kpc, with a mass ratio in the range ∼2 : 1. Dynamical estimates (Menanteau et al. 2012) suggest a high relative velocity for the merging infall, in the range ∼1500 km s⁻¹ to ∼2500 km s⁻¹, with the uncertainties in the allowed range 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 ∼ 600 kpc, and the SZ centroid is offset by a distance of ∼150 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). It is located in the SE cluster component and is offset by ∼70 kpc with respect the SE X-ray peak.

Moreover, the X-ray morphology of El Gordo presents several interesting features. The total X-ray luminosity is L_X ∼ 2 × 10⁴⁵ ergs⁻¹ in the 0.5 − 2 keV band (Menanteau et al. 2012); the highest X-ray emission is in the SE region and is 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 (Lindner et al. 2014).

This set of multifrequency observations led to the conclusion that El Gordo is a cluster at high redshift that is experiencing a major merger of two massive clusters that collided at a high velocity (≳2000 km s⁻¹). 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 the NW to the SE, inducing a wake in the ICM of the primary and the two-tail cometary structure seen in the X-ray images. The original gas core of the primary is absent because it was destroyed during the collision with the high-velocity compact core of the secondary. The secondary is thus 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 ΛCDM paradigm. For example, Asencio et al. (2021) argue that in a ΛCDM context the detection of a cluster such as 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 is instead 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 (unlike just after the pericenter passage); instead, 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 but appears to also be spatially offset from the mass centroid. This behavior has been observed in other clusters (Hallman & Markevitch 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). Simulations by 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 and the SZ centroid, but failed to replicate the observed twin-tailed X-ray morphology. This could be due to a number of factors, including the numerical scheme used in the simulations. Donnert (2014) carried out the simulations by employing a smoothed particle hydrodynamics (SPH; Price 2012) scheme, which, in its standard formulation, it is 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 that is higher (≳40°) 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, while 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 that is best able to reproduce many of the observed properties of El Gordo is a merging cluster (“model B”) with total mass ∼3 × 10¹⁵ M_⊙ and a high mass ratio (∼3.6). The initial merging configuration is that of an off-axis merger with impact parameter ∼800 kpc and a very high initial relative velocity (∼2500 km s⁻¹). The best match with measured offsets and X-ray images of El Gordo, as well as its temperature and luminosity, is obtained t ∼ 0.14 Gyr after it passes the pericenter and with a modest inclination angle (∼30°) 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 (∼5%). This in turn produces X-ray emission in the outer region of the merging cluster that is lower than observed. Finally, it is worth noting that the ZYL15 simulations can reproduce the total X-ray luminosity of El Gordo only for a primary cluster mass of ∼2.5 × 10¹⁵ M_⊙. 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 light of recent works that suggest substantially lower 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 on the order of ∼10¹⁵ M_⊙. This value is about a factor of 2 lower then previous SL studies (Zitrin et al. 2013). In a recent paper, Kim et al. (2021) present an improved WL analysis aimed at estimating the total mass of El Gordo and that of its components. The values they report for the NW and SE clusters are ∼9.9 × 10¹⁴ M_⊙ and ∼6.5 × 10¹⁴ M_⊙, respectively. These values are significantly lower (∼60%) than those dynamically inferred by Menanteau et al. (2012). They are also lower (∼30%) than the WL values reported by Jee et al. (2014).

These updated mass estimates were obtained independently by both SL (Diego et al. 2020) and WL studies (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 ΛCDM model (Kim et al. 2021; but see 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 of a merging scenario that incorporates the recently revised El Gordo cluster masses. The simulations were performed for a variety of merging initial conditions: we did this in order to find the final merging configuration that can best reproduce various observations of El Gordo. The initial condition setup mirrors that implemented by ZYL15; this was purposely chosen to ease the comparison of our results with previous findings.

The structure of the paper is as follows. 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 Sect. 3, with Sect. 3.1 showing results from off-axis merger simulations aimed at reproducing the observed twin-tailed X-ray morphology. In Sect. 3.2 we analyze the consistency of a returning scenario according to the results extracted from head-on cluster collisions, while Sect. 3.3 is dedicated to merger simulations that incorporate a BCG stellar component and to the impact of SIDM on some merging models. Section 4 summarizes our main conclusions. Throughout this work we use a concordance ΛCDM cosmology, with Ω_m = 0.3, Ω_Λ = 0.7, and Hubble constant H₀ = 70 ≡ 100 h km s⁻¹ Mpc⁻¹.

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, while 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 refer to the SPH scheme used here as integral SPH (ISPH).

For the sake of completeness, we introduce here 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 r_i, and velocity v_i that represent the fluid.

The SPH gas density ρ(r) at the particle position r_i

$$\begin{array}{c}\hfill \rho ({\mathit{r}}_{\mathit{i}})={\displaystyle \sum _j}{m}_{j}W(|{\mathit{r}}_{\mathit{ij}}|,h_i)\end{array}$$(1)

gives the particle gas density ρ(r_i)≡ρ_i. The summation should be over all the set of N particles, but because the adopted kernel W(|r_ij|,h_i) has compact support it reduces to only the N_nn neighboring particles j for which |r_i − r_j|≤2h_i. We use here the M₄ kernel (Price 2012). The smoothing length h_i is obtained by finding the root of

$$\begin{array}{c}\hfill \frac{4\pi {(2h_i)}^3{\rho }_i}{3}=N_{\mathit{nn}}{m}_i,\end{array}$$(2)

for which we set N_nn = 32 ± 3.

2.2. Initial conditions

To study the merging cluster El Gordo, we performed a series of N-body/hydrodynamical ISPH simulations. The merging system consists of a primary cluster with mass M₁ and a secondary with mass M₂, and the mass ratio is q = M₁/M₂.

The cluster mass is defined as the mass M₂₀₀ within the cluster radius r₂₀₀. This is the radius r_Δ within which the total mass is

$$\begin{array}{c}\hfill M_{\mathrm{\Delta }}=\frac{4\pi }{3}\mathrm{\Delta }{\rho }_c(z)r_{\mathrm{\Delta }}^3,\end{array}$$(3)

and the mean density is Δ times the cosmological critical density ρ_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 used for the two colliding clusters follows the same setting adopted by ZYL15. This was 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 performed a particle realization of two individual spherical halos, initially in hydrostatic equilibrium. Each halo consists of multiple mass components: DM, gas, and eventually star particles.

2.2.1. Dark matter

For the DM density profile we assumed an NFW profile

$$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}(r)=\frac{{\rho }_s}{r/r_s{(1+r/r_s)}^2},0\le r\le r_{200},\end{array}$$(4)

where ρ_s and r_s are the scaling parameters for density and radius, respectively. The scale radius r_s is defined by r_s = r₂₀₀/c₂₀₀, where c₂₀₀ is the concentration parameter. We set the value of c₂₀₀ using the c − M relation of Duffy et al. (2008):

$$\begin{array}{c}\hfill c_{200}=5.71a^{0.47}{\left(\frac{M_{200}}{2\times {10}^{12}{M}_{\odot }{h}^{-1}}\right)}^{0.084},\end{array}$$(5)

where a = 1/(1 + z). The DM cluster density profile is then specified once the mass M₂₀₀ and the concentration parameter c₂₀₀ are given.

At radii larger than r₂₀₀ the cumulative mass profile diverges as r → ∞. We therefore introduce for r >  r₂₀₀ an exponential cutoff up to a final radius r_max = νr₂₀₀:

$$\begin{array}{cc}& {\rho }_{\mathrm{DM}}(r)={\rho }_{\mathrm{DM}}(r_{200}){(r/r_{200})}^{\delta }exp{\displaystyle -\left(\frac{r-r_{200}}{r_{\mathrm{decay}}}\right)},\hfill \\ \hfill & r_{200}<r<r_{\max },\hfill \end{array}$$(6)

where r_decay = ηr₂₀₀ is the truncation radius, and the parameter δ is set by requiring the first derivative of the DM density profile to be continuous at r = r₂₀₀. For the runs presented here, we set the truncation parameters to the values (ν, η)=(2, 0.2); this choice was motivated due to 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 computed the cumulative DM mass M_DM(< r) within the radius r. We solved 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 ρ_m(r) of a given mass component is given by

$$\begin{array}{c}\hfill {\rho }_m(r)=4\pi {\displaystyle {\int }_0^{\mathrm{\Psi }(r)}}{f}_m(\mathcal{E})\sqrt{2[\mathrm{\Psi }(r)-\mathcal{E}]}d\mathcal{E},\end{array}$$(7)

where Ψ(r)= − Φ(r) is the relative gravitational potential of the system, ℰ = Ψ − v²/2 is the relative energy and f_m(ℰ) 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 Eq. (7) to obtain

$$\begin{array}{c}\hfill f_m(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}\left[{\displaystyle {\int }_0^{\mathcal{E}}}\frac{d^2{\rho }_m}{d{\mathrm{\Psi }}^2}\frac{d\mathrm{\Psi }}{\sqrt{(\mathcal{E}-\mathrm{\Psi })}}\right],\end{array}$$(8)

where the second-order derivative d²ρ_m/dΨ² can be evaluated analytically from ρ_m(r) and M(< r).

We constructed tables of f_m(ℰ) by evaluating numerically the integral (8) over a range of energies. We then used an acceptance–rejection method to obtain the particle speed $v=\sqrt{2[\mathrm{\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.

2.2.2. Baryonic matter

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

$$\begin{array}{c}\hfill {\rho }_{\mathrm{gas}}(r)=\frac{{\rho }_0}{(1+r/r_c)[1+{(r/r_c)}^2]},0\le r\le r_{200},\end{array}$$(9)

where r_c is the gas core radius and ρ₀ the central gas density. At radii larger than r₂₀₀ the density profile is continued by assuming that it follows that of the DM (Zhang et al. 2014):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{gas}}(r)={\rho }_{\mathrm{DM}}(r)\frac{{\rho }_{\mathrm{gas}}(r_{200})}{{\rho }_{\mathrm{DM}}(r_{200})},r_{200}<r<r_{\max }.\end{array}$$(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 considered the commonly employed non-isothermal β-model (Cavaliere & Fusco-Femiano 1978):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{gas}}(r)={\rho }_0{(1+\frac{r^2}{r_c^2})}^{-\frac{3}{2}\beta }.\end{array}$$(11)

For a given gas core radius, r_c, and slope parameter, β, to completely specify the gas density profile, it is necessary to determine the central density ρ₀. This is done by assigning the initial gas mass fraction f_g of the cluster at r = r₂₀₀ and then finding the corresponding root value ρ₀. The choice of the initial gas density profiles is a crucial issue in shaping the final X-ray morphology of the merging system. We provide a thorough discussion of the adopted profiles in Sect. 2.5.

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

$$\begin{array}{c}\hfill T(r)=\frac{\mu m_p}{k_B}\frac{G}{{\rho }_{\mathrm{gas}}(r)}{\displaystyle {\int }_r^{\infty }}\frac{{\rho }_{\mathrm{gas}}(u)}{u^2}{M}_{\mathrm{tot}}(<u)\mathrm{d}u,\end{array}$$(12)

where M_tot(< r) is the total mass within the radius r and we assume μ = 0.59. The gas thermal energy in the SPH equations is obtained by setting γ = 5/3 for the ideal equation of state.

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

2.2.3. Stellar component

For some of our merging runs, we also considered halos that contain a star matter component, in addition to the DM and gas. This is supposed to mimic the BCG mass distribution, which we initially placed at the center of the halo.

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

$$\begin{array}{cc}\hfill {\rho }_{\star }(r)& ={\rho }_{e}exp{\displaystyle \{-d_{n_s}[{(r/r_e)}^{1/n_s}-1]\}},\hfill \\ \hfill d_{n_s}& \sim 3n_s-1/3+0.0079/n_s,\mathrm{for}{n}_s\gtrsim 0.5,\hfill \end{array}$$(13)

where r_e is the effective radius and ρ_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):

$$\begin{array}{c}\hfill {log}_{10}{M}_{\star ,\mathrm{BCG}}\sim 0.38[{log}_{10}(M_{500})-14.5].\end{array}$$(14)

For M₅₀₀ ∼ 10¹⁵ M_⊙, which is appropriate for the range of halo masses considered here, the relation (14) gives M_⋆, BCG ∼ 2.3 × 10¹² M_⊙. We can now obtain the density ρ_e by solving for the ratio M_⋆, BCG/M₅₀₀ within r₅₀₀, once r_e is given.

If one adopts the relation of Shen et al. (2003, Eq. (33)), then for M_⋆, BCG ∼ 3 × 10¹² M_⊙ one obtains r_e ≃ 30 kpc. The smaller is the value of the effective radius, the higher must be the central density for a given M_⋆, BCG. However, the mass resolution of the simulations sets a lower limit for the minimum spatial scale that can be resolved, so 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 ∼10¹⁵ M_⊙ halo mass and M_⋆, BCG ∼ 3 × 10¹² M_⊙, we obtain ρ_⋆(r ∼ 10⁻²r₂₀₀)≃10⁴ρ_c ≃ 3.7 × 10⁶ M_⊙ kpc⁻³ for the central value of the stellar mass density. For a given initial star density distribution, we then determined positions and velocities of the star particles according to the procedure described in Sect. 2.2.1.

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

$$\begin{array}{c}\hfill m_d\simeq 8\times {10}^8{M}_{\odot }(M_{200}^{(2)}/2\times {10}^{14}{M}_{\odot })\end{array}$$(15)

and m_g = f_bm_d/(1 − f_b)) ≃ 0.16m_d, respectively. Here $M_{200}^{(2)}$ is the mass of the secondary and f_b ≃ 0.16 the cosmological mass gas fraction. With this mass resolution we have N_DM ≃ 3.4 × 10⁵ DM particles for an halo mass of M₂₀₀ ∼ 10¹⁵ M_⊙, and N_g ≃ 1.7 × 10⁵ gas particles for an halo gas mass of M_g ∼ 10¹⁴ M_⊙.

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)

$$\begin{array}{c}\hfill {\epsilon }_i=15.8·{(m_i/6.2\times {10}^8{M}_{\odot })}^{1/3}\mathrm{kpc}.\end{array}$$(16)

For 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 applied Eq. (15), but with a coefficient scaled down by a factor of 4. This is to prevent any discreteness effects on the star particle distribution during the simulations. With this choice of parameters the lower limit ${\epsilon }_{\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 Eq. (16). We note, however, that this lower limit is derived by requiring discreteness effects to be negligible over a Hubble time, whereas here the timescale of the collisions is ≲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, Δt_i, which is chosen from a power-of-two hierarchy Δt_m = Δt₀/2^m, where m = 0, 1… and Δt₀ = 1/100 Gyr is the largest allowed timestep. During the simulations, we recorded particle properties at simulation times t_n = nΔt₀, where n = 0, 1, … is the simulation step.

2.2.4. Initial merger kinematics

To construct the orbits of our merger simulations, we chose 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 {Xⁱⁿ, Vⁱⁿ} of the two clusters are set as in ZYL15. For a given mass ratio q = M₁/M₂ the initial separations are then X₁ = {d_ini/(1 + q),P/(1 + q),0} and X₂  =  {−d_iniq/(1 + q), − Pq/(1 + q),0} for the primary and secondary cluster, respectively. Here P is the impact parameter of the collision and $d_{\mathrm{ini}}=2(r_{200}^1+r_{200}^2)$. Similarly, for the initial velocities we set V₁ = { − V/(1 + q),0, 0} for the primary and V₂ = {Vq/(1 + q),0, 0} for the secondary, where V is the initial collision velocity.

We then shifted the particle positions and velocities of the two halos according to the vectors {Xⁱⁿ, Vⁱⁿ}; the merger dynamical evolution was fully determined by the merging parameters {M₁,  q,  P,  V}. In terms of these parameters the energy of the system can be expressed as

$$\begin{array}{c}\hfill E=\frac{1}{2}{M}_1{V}^2[\frac{1}{1+q}-\frac{G{M}_1}{d_{\mathrm{ini}}}\frac{1}{V^2}\frac{2}{q}],\end{array}$$(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₂₀₀.

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, SIDM 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 that represents a patch of DM phase-space that 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 restricted ourselves to the simplest case of isotropic and elastic scattering between DM particles; moreover, we further assumed a constant, velocity-independent DM cross-section σ_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), it can nonetheless 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 that we followed here was introduced by Vogelsberger et al. (2012), and its implementation requires the definition of a local DM density.

For each DM particle, i, we defined a DM density ρ_DM(r_i) according to the SPH formulation, and, as in Eq. (1), we thus defined a DM smoothing length $h_i^{\mathrm{DM}}$ such that

$$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}({\mathit{r}}_{\mathit{i}})={\displaystyle \sum _j}{m}_j^{\mathrm{DM}}W(|{\mathit{r}}_{\mathit{ij}}|,h_i^{\mathrm{DM}}),\end{array}$$(18)

where $m_i^{\mathrm{DM}}$ denotes the mass of the DM simulation particle i and the summation is over N_nn = 32 ± 3 DM neighboring particles. The DM smoothing length $h_i^{\mathrm{DM}}$ is evaluated whenever the particle i is active and its timestep is synchronized with the current simulation time. For a DM particle that 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 Δt_i is then determined as Vogelsberger et al. (2012)

$$\begin{array}{c}\hfill P_{\mathit{ij}}=m_i^{\mathrm{DM}}W(r_{\mathit{ij}},h_i^{\mathrm{DM}})\frac{{\sigma }_{\mathrm{DM}}}{m_X}{v}_{\mathit{ij}}\mathrm{\Delta }{t}_i,\end{array}$$(19)

where v_ij = |v_i − v_j| is the relative velocity between particles i and j, m_X is the physical mass of the DM particle and Δt_i the timestep of particle i.

The total probability of particle i scattering is P_i = ∑_jP_ij/2, where the factor 2 accounts for the other member of the scattering pair. A collision of particle i with one of its neighbors will occur if P_i ≤ 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 that satisfies $x\le {\displaystyle {\sum }_j^k}{P}_{\mathit{ij}}$.

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

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\mathit{u}}_{\mathit{i}}\hfill & =\mathit{V}+(v_{\mathit{ij}}/2)\mathit{e}\hfill \\ {\mathit{u}}_{\mathit{j}}\hfill & =\mathit{V}-(v_{\mathit{ij}}/2)\mathit{e},\hfill \end{array}\end{array}$$(20)

where V = (v_i + v_j)/2 is the center-of-mass velocity, and 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; it is described in Appendix A. We also tested the scattering implementation by showing results from an isolated SIDM halo.

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 σ_DM = 0 simulations. We justified this assumption by assuming the impact of initially cored DM density profiles to be subdominant in the evolution of post-pericenter gas structures (ZuHone et al. 2019).

2.4. Observables and projection analyses

For each of our merger simulations, we saved particle properties at various output times. This was done 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 could 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 adopted for the observer the same coordinate system $\{\widehat{x},\widehat{\mathit{y}},\widehat{z}\}$ introduced by the authors to construct their projected maps.

In such a system the line of sight is along the $\widehat{z}$ direction and the $\{\widehat{x},\widehat{\mathit{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 α, the second rotates the tilted frame by an angle i around the x′-axis. The angle i between the axes z and $\widehat{z}$ is thus the angle between the merging and sky plane. Our viewing direction is thus defined by the angles {α, i}, we refer to ZYL15 for more details on the setup of the rotation angles.

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

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_m(x,y)={\displaystyle {\int }_{\mathrm{los}}}[{\rho }_{\mathrm{gas}}(\mathit{x})+{\rho }_{\mathrm{DM}}(\mathit{x})]\mathrm{d}z,\end{array}$$(21)

where ρ_gas(x) and ρ_DM(x) are defined at the position x according to Eqs. (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¹ to their input X-ray maps. Here we followed a simpler approach and, as in Molnar & Broadhurst (2015), we evaluated the X-ray surface brightness by integrating the X-ray emissivity ε(ρ_g, T_g, Z, ν) along the line of sight and over the energy:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_X(x,y)& =\frac{1}{4\pi {(1+z)}^4}{\displaystyle {\int }_{\mathrm{los}}}\mathrm{d}z\hfill \\ \hfill & {\displaystyle \int \epsilon }({\rho }_g,T_g,Z,\nu )A_{\mathrm{eff}}(\nu )\mathrm{d}\nu ,\hfill \end{array}$$(22)

where T_g is the gas temperature, ν the frequency, Z the metal abundance of the gas, and A_eff(ν) the effective area of the telescope. For this, we used the four-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.3 Z_⊙; for reference, Z_⊙ ∼ 0.02 (Anders & Grevesse 1989) is the solar value.

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

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

$$\begin{array}{c}\hfill T_W=\frac{\int T(\mathit{x})W{d}^{3}x}{\int W{d}^{3}x},\end{array}$$(23)

where W is the weight function. We considered for W(x) the commonly employed gas density function W = ρ_g (mass-weighted temperatures, T_mw) and the weight function $W={\rho }_{\mathrm{g}}^2{T}^{-3/4}$, which defines the corresponding spectroscopic-like temperature T_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 line of sight.

The SZ surface brightness is calculated at frequency ν and including relativistic corrections (Itoh et al. 1998):

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\mathrm{SZ}}(x,y)& =\frac{{\sigma }_T{k}_B}{m_e{c}^2}{\displaystyle {\int }_{\mathrm{los}}}{n}_e{T}_g\hfill \\ \hfill & [g(\nu )+{\mathrm{\Sigma }}_{k=1}^{k=4}{Y}_k{\mathrm{\Theta }}^k]\mathrm{d}z,\hfill \end{array}$$(24)

where m_e, k_B, σ_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(ν)=coth(x_ν/2)−4 is the nonrelativistic frequency function, where x_ν = h_Pν/(k_BT_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 given in Itoh et al. (1998) and Θ ≡ k_BT_g/m_ec². As in ZYL15 we set ν = 150 GHz and smooth the SZ maps with a Gaussian kernel with width σ_SZ = 270 kpc (∼0.55′ at z = 0.87).

We evaluated 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 x_g is the sum

$$\begin{array}{c}\hfill f({\mathit{x}}_g)={\displaystyle \sum _i}{f}_i\frac{m_i}{{\rho }_i}W(|{\mathit{x}}_g-{\mathit{x}}_i|,h_i)\end{array}$$(25)

over the subset of SPH particles that satisfy |x_g − x_i|≤2h_i. To perform the integrals along the line of sight, we exploited the property of the M_n splines (Price 2012), which approximate a Gaussian in the limit n → ∞. For the M₄ kernel used here, the integrals along the line-of-sight z were then performed analytically, and we calculated projected quantities at the grid point x_g as

$$\begin{array}{c}\hfill f({\mathit{x}}_g)={\displaystyle \sum _i}{f}_i\frac{m_i}{{\rho }_i}{W}_{2\mathrm{D}}(x_{\mathit{gi}}^2,{\sigma }_i^G),\end{array}$$(26)

where $x_{\mathit{gi}}^2={(x_g-x_i)}^2+{({\mathit{y}}_g-{\mathit{y}}_i)}^2$, ${\sigma }_i^G\simeq h_i/\sqrt{3}$ (Dehnen & Aly 2012) and W_2D is the 2D Gaussian,

$$\begin{array}{c}\hfill W_{2\mathrm{D}}=\frac{1}{{(\sqrt{2\pi }{\sigma }_i^G)}^2}exp{[-x_{\mathit{gi}}^2/2{\sigma }_i^G]}^2.\end{array}$$(27)

To locate the centroid positions of the various maps that we produced to study the El Gordo cluster, we applied the shrinking circle method to the simulation particles. 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 gas, DM, or star particles that 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 x_c = ∑_iw_ix_i/∑_iw_i, where the weights w_i are chosen according to the centroid under consideration.

To locate the centroid position, we then started from an initially large radius R⁽⁰⁾(∼5 Mpc) that contains all of the particle subset, and we calculated the weighted mean position, ${\mathit{x}}_c^{(0)}$. We then proceeded to iteratively calculate at each step k = 1, 2, …, M the centroid position ${\mathit{x}}_c^{(k)}$, based on the subset particles that are within a sphere of radius R^(k) = fR^{(k − 1)} <  R^{(k − 1)}, located at ${\mathit{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 ${\mathit{x}}_c^{(M)}$. We found that the results are robust for N_M ∼ 100, but the choice of the shrinking factor, f, is critical: it must be chosen sufficiently large (f ≳ 0.8) to prevent the algorithm from failing when there are multiple peaks.

The method is easily applied in 2D to the projected particle positions along the viewing direction. Therefore, for a given mock image, we calculated 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 that 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 that 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 ∼ 700 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: {α, i}={−90° ,30° }.

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 Fig. 14 of Menanteau et al. (2012) that the map of the deprojected electron density exhibits a peak value of n_e ≃ 5 × 10⁻² cm⁻³ 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, that 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 constructed our ensemble of off-axis mergers by considering different values of the primary mass $M_{200}^{(1)}$ and modifying the mass ratio q = M₁/M₂ such that the mass of the secondary stays close to the estimated value $M_{200}^{(2)}\sim 6.5\times {10}^{14}{M}_{\odot }$. From Eq. (5) one obtains ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2.682$ for the halo concentration parameter of the secondary. In accordance with recent mass estimates (Diego et al. 2020; Kim et al. 2021), we examined mergers with the following masses for the primary: $M_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$ and $M_{200}^{(1)}={10}^{15}{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 that we analyze, together with the corresponding merging parameters $\{M_{200}^{(1)},q,P,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 Sect. 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_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$, panel f in their Fig. 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_{200}^{(1)}\simeq 2.5\times {10}^{15}{M}_{\odot }$) is required in order to obtain an X-ray luminosity as high as that observed in the El Gordo cluster.

We argue that a viable solution for obtaining consistency between the recently measured cluster masses and observed luminosities is to consider the possibility of an initial gas density profile for the primary that is different from the Burkert one (9) previously adopted. Specifically, we employed the well-known β-model (11) to describe the initial radial gas density profile of the primary. For a given central density, ρ₀, the radial behavior of the gas density can be modified by varying the gas scale radius, r_c, and the slope parameter, β, with the aim of constructing the initial conditions for a merging system that can satisfy the observational constraints previously described.

We further simplified the parameter space of our merger models by setting β = 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_{200}^{(1)},q,P,V\}$, this choice leaves us with the gas core radius, r_c, of the primary as the only parameter needing to be set up for the initial conditions of our merger models.

We then constructed 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_{200}^{(1)},q,V,P\}$, we set up the initial conditions of our models by considering a range of values for the core radius, r_c. Table 2 lists the merger models that we considered; the ID of each model is 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 ζ = r_s/r_c. Hereafter, the initial gas core radius of the primary will be denoted as $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$.

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

3. Results

This section is dedicated to the presentation of our main results, which we divide into three subsections. In Sect. 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 Sect. 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 Sect. 3.3 some of the merging runs of Sect. 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 are 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 Sect. 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 Fig. 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 {α, i} of the merging plane, as in ZYL15. The bottom panels (model A_1) can then be compared directly with those in Fig. 1 of ZYL15, while the top panels (model B_1) correspond to those of their Fig. 3. Each panel refers to the post-pericenter epoch when the projected distance between the mass centroids is d_DM ∼ 700 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 Fig. 1, which exhibit the same morphological features seen in Figs. 1 and 3 of ZYL15. In particular, for model A_1 the bottom-left ({α, i}={−50° ,0° }) and bottom-right ({α, i}={0° ,0° }) panels show the same X-ray structures seen in Figs. 1a and c of ZYL15: a strongly asymmetric tail (Fig. 1a) and a small asymmetric wake-like feature (Fig. 1c). The bottom-middle panel ({α, i}={−50° ,75° }) of Fig. 1 corresponds to the fiducial model A of ZYL15 (Fig. 1b). The X-ray morphology is more symmetric here and very similar to that seen in Fig. 1b, with a bullet-like shape.

Similarly, the top panels of Fig. 1 can be compared with those displayed in Fig. 3 of ZYL15 for their model B. The top-left ({α, i}={−90° ,0° }) and top-right ({α, i}={−90° ,60° }) 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 Figs. 3a and c of ZYL15. The top mid panel of Fig. 1 corresponds to Fig. 3b of ZYL15 and to their fiducial model B. A comparison with this model is particularly significant, since this is the model that best matches the observations. In this panel ({α, i}={−90° ,30° }) the X-ray morphology is clearly characterized by a well-defined twin-tailed structure, as seen in observations and in their Fig. 3b.

The consistency between our model B_1, with viewing angles {α, i}={−90° ,30° }, 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_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_X we use Eq. (23).

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 ≃ 2.86 × 10⁴⁵ ergs⁻¹ and T_X ≃ 12 keV, respectively. These values can be directly compared with the corresponding ones reported in Table 2 of ZYL15: L_X ≃ 2.05 × 10⁴⁵ ergs⁻¹ and T_X ≃ 15 keV. From the same table one also has δT ≃ −1.1 × 10^−3 °K for the central SZ decrement, while from Table 4δT ≃ −9 × 10^−4 °K. 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.

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 accordance with the corresponding fiducial model B of ZYL15. We therefore adopted model B_1, which we 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, 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 Fig. 1. The main collision parameters that 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 that is produced after the collision thus depends on the adopted set of merging parameters {M₁,  q,  P,  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 Sect. 3.1.2). 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. Figs. 5a and b of ZYL15).

We thus deduce that for a given set of merging parameters {M₁,  q,  P,  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_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$: Models Bf, Bg, and Bh

We next used our 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 considered two families of off-axis merger models. The first family has the mass of the primary set to $M_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$ (models Bf, Bg and Bh), while $M_{200}^{(1)}={10}^{15}{M}_{\odot }$ for the second family (models Bk, Bl and Bj). This set of models will be presented in Sect. 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_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$. This is about ∼50% higher than the mass estimate (∼9.9 × 10¹⁴ M_⊙) 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 that progressively deviate from those of the reference model B. The first model that we consider is model Bf, with P = 600 kpc and V = 2500 km s⁻¹, followed by models Bg (V = 2000 km s⁻¹) and Bh (V = 1500 km s⁻¹).

Model Bf is nothing else than the merger model presented by ZYL15 in their Fig. 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₁,  q,  P,  V}, can reproduce the observational features of the El Gordo cluster.

As introduced in Sect. 2.5, we adopt here the simulation strategy of performing a set of merging simulations with initial conditions that differ in the choice of the primary’s gas scale radius $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$, 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 Fig. 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. The current epoch is defined when the projected distance d_DM between the mass centroids is approximately d_DM ∼ 700 kpc, and as for the fiducial model B of Fig. 1 the viewing angles are {α, i}={−90° ,30° }. 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 r_s ≃ 0.7 Mpc, so the $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$ ranges from $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 145\mathrm{k}\mathrm{p}\mathrm{c}$ (ζ = 4.81: Bf_rc14) up to $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 290\mathrm{k}\mathrm{p}\mathrm{c}$ (ζ = 2.4: Bf_rc29). From the maps of Fig. 2 it appears then that a change of $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$, 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.

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 Fig. 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^{\mathrm{p}\mathrm{r}\mathrm{m}}$ 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 Fig. 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.

It should also be stressed that such a dependence of the amount of stripped gas on $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$ takes place because the gas mass fraction is held fixed between the models. This in turn implies that the central gas density of the primary is anticorrelated with $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$. For a central gas density kept constant, one expects the dependence of the amount of stripped material on $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$ to be reversed.

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^{\mathrm{p}\mathrm{r}\mathrm{m}}$. Model Bf_rc20 (top-right panel of Fig. 2) is the model whose X-ray image is best able to match the requested morphology: for this model $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 200\mathrm{k}\mathrm{p}\mathrm{c}$ and its twin-tailed morphology well reproduces that of fiducial model B in Fig. 1. Hereafter, we denote the fiducial model for the Bf family of merger models as model Bf_rc20.

The same reasoning can be applied to the Bg models, for which the bottom panels of Fig. 2 show the mock X-ray maps extracted from three runs with different values of the gas core radius of the primary, $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$. For these models the initial collision velocity is lower (V = 2000 km s⁻¹) than in the Bf models, and model Bg_rc20 $(r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 200\mathrm{k}\mathrm{p}\mathrm{c})$ is now the fiducial model that produces the optimal X-ray morphology. We 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 3 shows for the Bh models X-ray maps extracted from three merger simulations. For these models V = 1500 km s⁻¹ and it is interesting to note that all of the considered models, from $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 290\mathrm{k}\mathrm{p}\mathrm{c}$ (ζ = 2.4: Bh_rc29) down to $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 146\mathrm{k}\mathrm{p}\mathrm{c}$ (ζ = 4.81: Bh_rc14), produce strongly asymmetric X-ray structures that 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_{200}^{(1)}\simeq 1.6\times {10}^{15}{M}_{\odot }$, require collision velocities at least as high as V ≳ 2000 km s⁻¹ 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_{200}^{(1)}={10}^{15}{M}_{\odot }$: Models Bk, Bl, and Bj

For the family of merger models with mass of the primary $M_{200}^{(1)}={10}^{15}{M}_{\odot }$, we have model Bk (P = 600 kpc and V = 1500 km s⁻¹) and Bl (V = 2000 km s⁻¹), 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 Fig. 4 demonstrate that the observed twin-tailed morphology cannot be reproduced by merger models with initial collision velocity as low as V = 1500 km s⁻¹ (Bk runs), regardless of the chosen value of the primary’s gas core radius. This behavior is similar to that of model Bh (Fig. 3), but now the mass of the primary is lower.

For Bl models, the collision velocity is increased to V = 2000 km s⁻¹. The bottom panels of Fig. 4 show that this class of merger model can produce an X-ray image that compares well visually with those previously used to define the list of fiducial models. Model Bl_rc24 (bottom-left panel of Fig. 4) is the merger model that is best able to match the observed X-ray morphology: for this model r_s ≃ 0.58 Mpc and $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 240\mathrm{k}\mathrm{p}\mathrm{c}$ 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 essentially considered only merger models with initial impact parameter P = 600 kpc; the exception is fiducial model B (P = 800 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 kpc down to P = 400 kpc. For each set of collision parameters {V,  P} we performed merging simulations by considering different values of $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}$. 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_{200}^{(1)}={10}^{15}{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 kpc are unable to generate the requested X-ray morphology. The X-ray structure is strongly asymmetric, with only one tail present if V = 2000 km s⁻¹ (Bjb_rc14 run) and approaches that of a head-on collision when V = 1500 km s⁻¹ (Bja_rc14 run). The asymmetric one-tail structure is even more pronounced if one considers model Bjc_rc12 for which P = 500 kpc.

3.1.4. Comparison of the results with previous findings and observations

In order to better constrain the fiducial merger models that 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_X given in Table 4 are obtained by using the estimator given by Eq. (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_X ≃ 12 keV) is smaller than that reported by ZYL15 for their fiducial model B (T_X ≃ 15 keV) by a factor of ∼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_X ≃ 13.5 keV) fare better than model B to approximate the observational estimate (T_X ≃ 14.5 keV). Model Bl_rc24 has the lowest temperature (T_X ≃ 11 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 dependence 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 identified these regions by visual inspection, with their spatial location and orientation 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 Δ across the wake; these can be compared with those depicted in Fig. 8 of ZYL15 or Fig. 4 of Molnar & Broadhurst (2015), after appropriate rescaling. All of the profiles agree quite well with each other, they share a minimum at Δ ∼ 0 arcsec and their two peaks are separated by about ∼50 arcsec. The profiles are approximately symmetric around Δ = 0, but exhibit a weaker emission in the wings, with respect the ones presented by ZYL15 and Molnar & Broadhurst (2015). The origin of such a discrepancy is not clear, we suggest that could be due to possible differences in the procedures adopted to define the extraction region used to calculate the spectra.

The separation between the two maxima and their emission level (∼1.5 × 10⁻⁸ counts s⁻¹ cm⁻² arcsec⁻²) are in accordance 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 of ∼2 lower. This discrepancy is less pronounced in the spectra of Molnar & Broadhurst (2015), for which the emission level of the minimum is about ∼70% of the maxima.

The SZ-X-ray offset (d_X − SZ ∼ 500 kpc) is substantial in accordance 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 (∼600 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 ≃ − 915 μ^°K, in line with observational estimates and in accordance with the corresponding value reported by ZYL15 for their model B (−1130 μ^°K). Models Bf_rc20 and Bg_rc20 have similar decrements, while for model Bl_rc24 the decrement is much smaller (≃ − 660 μ^°K) 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 that 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\text{km}{\text{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. In Sect. 3.3 we investigate the observational consequences when this assumption is no longer valid.

For the four fiducial merger models we show in Fig. 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 Fig. 9 of Menanteau et al. (2012). The figure also reports the difference $V_r^{\mathrm{DM}}$ 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 that, in the observer plane, are within a distance of 400 kpc from its mass centroid.

The velocities are in the range from $V_r^{\mathrm{DM}}\sim 1020\text{km}{\text{s}}^{-1}$ (model B) to $V_r^{\mathrm{DM}}\sim 850\text{km}{\text{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.

The two velocity histograms of fiducial model B can be directly compared with the corresponding ones shown in Fig. 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_r^{\mathrm{DM}}\sim 910\text{km}{\text{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 2000 km s⁻¹ ≲ V ≲ 2500 km s⁻¹ and 600 kpc ≲ P ≲ 800 kpc. The mass of the primary is allowed to vary in the range ${10}^{15}{M}_{\odot }\lesssim M_{200}^{(1)}\lesssim 2.5\times {10}^{15}{M}_{\odot }$, with the mass ratio q adjusted to satisfy $M_{200}^{(2)}\simeq 6.5\times {10}^{14}{M}_{\odot }$. The lower part of the mass interval is favored by the latest mass estimates.

The fiducial merger models that 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_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$, show an overall better agreement with data.

3.2. A returning scenario?

One of the most challenging features to reproduce 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, Fig. 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 considered head-on collisions and accordingly set P = 0; the procedures for setting up the initial conditions are those outlined in Sect. 2.

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

For a given value of q we considered three different values of the collision velocity: V = 1000,  2000, and 3000 km s⁻¹; our set of simulations thus 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 → 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 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 highly energetic (V = 3000 km s⁻¹), almost head-on collision. As discussed in the previous section, their fiducial model A corresponds to the bottom-middle panel of Fig. 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 Fig. 1 (Fig. 1c of ZYL15), which the authors interpret as a consequence of a nonzero impact parameter.

In this 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 ≃ 600 kpc, while we found here d_DM ≃ 720 kpc. We attribute this discrepancy to the very large value of the relative velocity (V_rel ≃ 4700 km s⁻¹) 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 Δt₀ = 1/100 Gyr. The current epoch is identified at run time when d_DM ≥ 700 kpc, and the snapshot of Fig. 1 (bottom right) corresponds to the simulation step n = 164. The previous output is saved at n = 162, when d_DM ≃ 630 kpc.

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 Fig. 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 ≃ 700 kpc after the pericenter passage; here t = 0 is when the simulations start.

From Fig. 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 R02V10, but its initial conditions are a bit peculiar (see below). 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 accordance with previous findings (Ricker & Sarazin 2001), while for higher mass ratios the post-pericenter shape of the ICM depends on the initial velocity V.

The only merging simulations that 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 the 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, Fig. 9 shows X-ray maps extracted from the same runs as in Fig. 8, but now in a post-apocenter phase. The criterion to identify the observer epoch is still that of having d_DM ≃ 700 kpc, but this when the secondary has reached the apocenter and is now returning toward the primary.

The time difference Δt_orb between the epochs shown in Fig. 9 and those in the corresponding panels of Fig. 8 gives, for a specific model, the time spent by the secondary between the outgoing phase and the returning one. The time Δt_orb can assume a wide range of values: from Δt_orb ∼ 0.02 (R02V10) up to Δt_orb ∼ 4 Gyr (R03V30).

The time Δ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 R02V10 has a very small value of Δt_orb because its initial velocity (V = 2000 km s⁻¹) is such that when the secondary reaches the orbit’s apocenter it is also when d_DM ≃ 700 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 had to choose from Fig. 9 the merging runs for which the value of the mass ratio is closer to the value they adopted (q ∼ 1.7). These conditions are satisfied by the runs R01V10, R01V20, and R01V30, for which Δt_orb ∼ 0.8, 0.4, and ∼1 Gyr, respectively. R01V30 is, however, ruled out by the timescale Δt_orb = 0.91 − 0.46 ∼ 0.5 Gyr favored by Ng et al. (2015). Their Monte Carlo simulations include neither dynamical friction nor tidal forces, so a comparison between their favored value of Δ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 ≳ 5). We then use these findings to argue that the constraint on Δt_orb of Ng et al. (2015) can be used to rule out model R01V30.

This leaves us with models R01V20 and R01V10. We 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 is approximately ∼100 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 Fig. 1c are short-lived. In particular, for model R01V20 the X-ray map of Fig. 9 can be compared with the maps presented by Donnert (2014, Fig. 5) in his simulation study of El Gordo. To properly compare the two merger models, it is necessary to identify among the various panels of Fig. 5 the post-apocenter epoch when the secondary is falling onto the primary, and the relative distance between the two clusters is approximately of ∼500 − 700 kpc. It is easily seen that this condition is satisfied at t ∼ 4.5 Gyr, and a comparison between the two panels reveals the presence of very similar X-ray structures.

The results presented here are valid for head-on collisions; nonetheless, it should be emphasized that off-axis mergers are also clearly inconsistent with a returning scenario. This can be clearly seen by realizing that for these merger models the time difference Δt_orb between the pre- and post-apocenter epochs, at which must be satisfied the constraint d_DM ≃ 700 kpc, are significantly higher than for mergers with orbits having zero angular momentum. For the fiducial models of Sect. 3.1 it is found Δt_orb ∼ 1 − 2 Gyr, whereas we have seen that the post-collision X-ray structures are fast evolving and their lifetime lies in the range ∼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 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 (∼100 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, nonzero offsets are also measured (Kim et al. 2021) between the BCG to DM centroid (∼60 kpc) and the BCG to X-ray peak (∼70 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 (∼10 000) Zitrin et al. (2012) found a relatively small mean offset (∼15 kpc) between the BCG and DM cluster centroids. More recently, De Propris 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 centers), we re-simulated two of the fiducial models outlined in Sect. 3.1. We chose fiducial models Bf_rc20 and Bl_rc24, whose simulation results produced the best fitting to the data overall and which have a mass of the primary that brackets the observational range.

Our new models are designated DBf_rc20 and DBl_rc24. The initial condition setup is the same as that of their respective model in Sect. 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 that are supposed to represent the BCGs is described in Sect. 2.2.3. The stellar masses of the primary and secondary cluster that we obtain from Eq. (14) are $M_{\star }^{(1)}=2.2\times {10}^{12}{M}_{\odot }$ and $M_{\star }^{(2)}=1.6\times {10}^{12}{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 Figs. 2 and 4, respectively.

About the dynamics of the BCGs two important results emerge from Fig. 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 Fig. 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 return to the likelihood of such a scenario for the merging cluster El Gordo in the conclusions.

An alternative scenario that could explain the observed spatial offsets between the different mass components involves SIDM that 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. 2017a, 2019; ZuHone et al. 2019; Fischer et al. 2022).

We investigate here the possible outcomes of SIDM cluster mergers, using the implementation of the SIDM modeling described in Sect. 2.3. We performed SIDM merging simulations for each of the two merger models of Fig. 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 SIDM cross-section assuming the following values: σ_DM/m_X = 1,  2, and 5 cm² gr⁻¹, 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 Fig. 11; the top panels are for model XDBf_rc20 and the bottom panels for model XDBl_rc24. We also considered SIDM simulations with σ_DM/m_X = 0.1 cm² gr⁻¹, 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 Fig. 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 Fig. 11 from the σ_DM/m_X = 5 cm² gr⁻¹ run of model XDBf_rc20 appear much rounder than those exhibited in Fig. 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 Fig. 11, this effect is present in all of the runs, but increases as σ_DM/m_X gets higher and when the mass of the primary gets lower. In fact, the bottom-right panel of Fig. 11 shows that for σ_DM/m_X = 5 cm² gr⁻¹ 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 (∼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 Sect. 3.2).

This effect can be clearly seen in the panels of Fig. 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 σ_DM/m_X = 1 cm² gr⁻¹ up to σ_DM/m_X = 5 cm² gr⁻¹, 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 Fig. 10, we also report in each panel the simulation step n. From this one can recover the simulation time t_n = nΔt₀, at which the maps were extracted.

We attribute this increasing dependence of the elapsed time, t, on the DM scattering cross-section, σ_DM/m_X, to the mutual transfer of energy that 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 σ_DM/m_X ∼ 5 cm² gr⁻¹. As a consequence of these dissipative processes, the post-pericenter cluster bulk velocities are now strongly reduced with respect the standard CDM model. As shown below, this in turn will also impact on the BCG bulk velocities.

Another important feature that 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 Fig. 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 Fig. 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 Fig. 11, there are several noteworthy features that 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 Fig. 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 Fig. 11 the relative separation between the two peaks initially decreases as σ_DM/m_X increases, and becomes positive for σ_DM/m_X = 5 cm² gr⁻¹. For these runs the X-ray centroid is no longer lagging behind the mass centroid but is now clearly ahead of it.

For the same value of σ_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 = 2000 km s⁻¹) in the former than in the latter (V = 2500 km s⁻¹). Leaving aside the σ_DM/m_X = 1 cm² gr⁻¹ runs, for which this dependence is too weak to be detected, the two merging simulations with σ_DM/m_X = 2 cm² gr⁻¹ (middle panels of Fig. 11) effectively reproduce this behavior.

However, for the merging runs with σ_DM/m_X = 5 cm² gr⁻¹ 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 ∼ 0.45 Gyr) than in model XDBf_rc20 (t ∼ 0.24 Gyr).

The same time difference is responsible of another noteworthy feature of model XDBl_rc24 with σ_DM/m_X = 5 cm² gr⁻¹. From the bottom-right panel of Fig. 11 it can be seen that this run exhibits negative BCG-DM offsets, that is, the DM peaks are farther 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. 2017a; Fischer 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 Fig. 11 for model XDBl_rc24. These behaviors are in line with what is expected in a merging SIDM scenario, due to the general nonlinear dependence on σ_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 σ_DM/m_X = 5 cm² gr⁻¹ is the only merging simulation that is able produce relative separations within the observational range (∼50 − 100 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 are in the range σ_DM/m_X ∼ 4 − 5 cm² gr⁻¹. 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 ∼ 700 kpc.

This is made clear by the middle panels of Fig. 11, where for σ_DM/m_X = 2 cm² gr⁻¹ both models exhibit a relative separation between the different centroids of about ∼50 kpc, at approximately the same epoch t. However, for σ_DM/m_X = 5 cm² gr⁻¹ 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 σ_DM/m_X → 5 cm² gr⁻¹ 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. We 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_{200}^{(1)}\sim 1.6\times {10}^{15}{M}_{\odot }$.

Motivated by these results, we now turn our attention to the SIDM merging simulation of model XDBf_rc20 with σ_DM/m_X = 5 cm² gr⁻¹, which from the top-right panel of Fig. 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 Fig. 1), we estimate for the SE cluster of this model a peak value of the projected DM mass surface density of about Σ_DM ∼ 0.3 gr cm⁻². The cluster’s scattering depth is then τ = Σ_DMσ_DM/m_X ∼ 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 τ.

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 Fig. 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 σ_DM/m_X = 5 cm² gr⁻¹.

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 σ_DM/m_X = 5 cm² gr⁻¹ 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_{200}^{(1)},q,V,P\}=\{1.6\times {10}^{15}{M}_{\odot },2.32,2500\text{km}{\text{s}}^{-1},600\text{kpc}\}$.

For this set of initial collision parameters we performed SIDM merger simulations for runs with both σ_DM/m_X = 5 cm² gr⁻¹ and σ_DM/m_X = 4 cm² gr⁻¹. 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 are required than for standard CDM runs.

For the σ_DM/m_X = 5 cm² gr⁻¹ SIDM merging run it was found after some testing that the best results are obtained for $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 400\mathrm{k}\mathrm{p}\mathrm{c}$ and cluster gas fractions (f_g1, f_g2)=(0.12, 0.12). We label this model as XDBf_sa. Similarly, for the σ_DM/m_X = 4 cm² gr⁻¹ simulation it was found that $r_c^{\mathrm{p}\mathrm{r}\mathrm{m}}\sim 300\mathrm{k}\mathrm{p}\mathrm{c}$ and (f_g1, f_g2)=(0.12, 0.14) are the optimal settings to match the observed X-ray morphology. We label this model as XDBf_sb. We report in Table 7 the most relevant initial merger parameters of the two models; in both cases the initial gas density profile of the secondary is the same as in the standard CDM runs. Figure 12 shows the mock X-ray maps extracted from these two SIDM merger simulations.

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

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

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

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

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

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

As already outlined, the most striking feature is the position for the SE cluster of the X-ray peak, which is now leading the DM mass centroid. Moreover, the relative distance d_X-DM between the two centroids varies from d_X-DM ∼ 60 kpc (σ_DM/m_X = 4 cm² gr⁻¹) up to d_X-DM ∼ 120 kpc (σ_DM/m_X = 5 cm² gr⁻¹). The agreement of this range of values with the measured offset $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}\sim 100\text{kpc}$ can be considered significant, given the observational uncertainties (see below).

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

Another important result that follows from an SIDM scenario in major cluster mergers is the expected presence of BCG-DM offsets, with BCG centroids that exhibit miscentred positions with respect to the centers of their parent DM halos. Figure 12 shows that according to value of σ_DM/m_X, the BCG-DM offsets for the SE cluster lie in the range d_BCG-DM ∼ 60 − 120 kpc, which agrees fairly well with the observed offset of about ∼60 kpc. We note that for the NW cluster, the position of the simulated BCG strictly coincides with that of the observed peak galaxy number density.

In SIDM mergers, positive BCG-DM offsets are expected after the first pericenter passage. This effect can be interpreted as a result of the deceleration of the DM component, this in turn is due to the exchange of energy between the DM halos that occurs during the cluster collision. As mentioned in the discussion of the results from the σ_DM/m_X = 5 cm² gr⁻¹ SIDM simulations of Fig. 11, positive BCG-DM offsets are a transient phenomena. As a BCG starts to climb out of the cluster potential after the pericenter passage, it will experience a gravitational pull (Harvey et al. 2014) from the DM mass component that will reduce its forward momentum as well as the size of the offset.

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

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

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

Moreover, the X-ray emission in the outer regions behind the secondary is significantly less than in model Bf_rc20. This can be clearly seen by the location of the third contour level (2.9 × 10⁻¹ counts arcsec⁻² from the inside), which in the X-ray map of model Bf_rc20 (top-right panel of Fig. 2) has a minimum along the line joining the two DM centroids at about ∼600 kpc from the X-ray peak of the secondary. For the merger simulation of model XDBf_sb this minimum is much less pronounced and at about ∼400 kpc from the centroid of the X-ray emission.

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

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

Finally, one of the most significant prediction of DM collisionality is the flattening of the DM core density. It is therefore important to assess in SIDM mergers the consistency of the predicted mass distribution with observed mass profiles. The scattering rate of DM particles is higher in high-density regions, and due to the collisions between DM particles the net effect is a transfer of energy toward the inner regions. This effect is evident from the left panel of Fig. A.1, in which for different values of σ_DM/m_X are shown at t = 1 Gyr the density profiles of an isolated SIDM halo. As predicted, the flattening of density at inner radii is higher as σ_DM/m_X increases.

The development of flatter DM density profiles is however expected to be significantly increased in SIDM cluster mergers. This is because the scattering probability (19) depends on the relative velocity between DM particles. At the pericenter the relative velocity between the two clusters can be as high as v_rel ∼ 4000 km s⁻¹, and this in turn implies much higher scattering rates than in the case of an isolated halo.

To explicitely verify this expectation, Fig. 13 shows for our model XDB_sb the DM radial density profiles of each cluster. For the sake of comparison these are depicted at the beginning of the simulation and after the core passage, at the observer epoch t = 0.24 Gyr. Additionally, at the same epochs we also show the density profiles of the BCG stellar component for the corresponding cluster. The profiles of the NW and SE cluster are plotted separately in the left and right panel, respectively.

As expected, after the core passage the inner density profile of the two DM halos exhibits a significant flattening toward the cluster center. Accordingly, for both of the DM density profiles we now find a core radius of approximately ∼200 kpc. It is interesting to make a comparison against observational estimates of the concentration parameters, as obtained by fitting the density profiles of each DM halo using an NFW model. We refer to Duffy et al. (2008) for a description of the adopted fitting procedure. From the best fit of each halo we obtain for the concentration parameters and their formal statistical errors ${\mathrm{c}}_{200}^{\mathrm{N}\mathrm{W}}=1.99\pm 0.026$ and ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=1.57\pm 0.024$ for the NW and SE cluster, respectively.

From their WL study on El Gordo cluster, in their Table 2 Kim et al. (2021) report for the two clusters ${\text{c}}_{200}^{\text{NW}}(K)=2.{54}_{-0.41}^{+0.97}$ and ${\text{c}}_{200}^{\text{SE}}(K)=3.2_{-0.74}^{+1.17}$. For the NW cluster the observational estimate ${\mathrm{c}}_{200}^{\mathrm{N}\mathrm{W}}(K)$ of the concentration parameter is marginally consistent at 1σ level with the best fit value ${\mathrm{c}}_{200}^{\mathrm{N}\mathrm{W}}$, as extracted from the SIDM merger simulation XDBf_sb. On the contrary, for the SE cluster the estimate ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}(K)\sim 3$ of Kim et al. (2021) is significantly higher than the best fit parameter ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}\sim 1.6$ predicted by our SIDM merger simulation.

The reason for such a large difference is in the adopted initial conditions of our merger simulations. For the SE cluster the initial value c₂₀₀ = 2.682 of the concentration parameter is predicted by the Duffy et al. (2008) relationship (Eq. (5)), and it was used in the merger simulations performed here to set up the initial conditions of the SE DM halo.

For the SE cluster this then implies that our initial setting of c₂₀₀ = 2.682 is much lower and marginally consistent, within the quoted uncertainties, with the observational estimate ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}(K)\sim 3$ of the DM halo concentration parameter. Therefore, it is not surprising that we find such a discrepancy between the value of ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}(K)$ and the best fit value ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}\sim 1.6$ extracted from the SIDM merger simulation XDBf_sb.

However, the presence of such a discrepancy should not be considered particularly significant. As pointed out by Kim et al. (2021), estimates of the concentration parameters are poorly constrained due to their low sensitivity to the lensing signal; moreover, the M − c relation (Eq. (5)) has been derived by averages from a simulation ensemble of DM halos, so its application to a post-collision system such as El Gordo should be considered with caution.

Additionally, in each panel of Fig. 13 we also show the BCG density profile of the corresponding cluster. The radial profile of each BCG is plotted at the present epoch (solid blue lines) and at t = 0 (dashed blue line). We can see that for both of the halos the two profiles exhibit very small differences, thus showing that the internal star distribution of the BCGs remains relatively unperturbed during the collision.

It is important to note that the profiles are centered on the mass centroids of the BCGs when they are calculated, so the location of their origin varies with time and a meaningful comparison between the BCG and DM density profiles can only be performed at t = 0. Using a sample of seven, massive, relaxed galaxy clusters Newman et al. (2013) showed that in cluster central regions (≲10 kpc) the star mass component becomes dominant.

There would therefore be a disagreement between these findings and the behavior at t = 0 of the radial density profiles of the two mass components. For each halo, Fig. 13 shows that at r ∼ 10⁻²r₂₀₀ ∼ 15 kpc the star density ρ_⋆ is lower than ρ_DM by a factor of ∼4. This discrepancy is partly due to a resolution effect, the effective radius of Sect. 2.2.3 being set to r_e = 60 kpc. Moreover, the profiles of Fig. 13 are presented in the chosen radial range of 0.01 ≲ r/r₂₀₀ ≲ 1, at the spatial resolution limit of r = 10^−2.8r₂₀₀ ∼ 3 kpc the star density is found to be higher than the DM density: ρ_⋆ ∼ 8 × 10⁵ρ_c ∼ 4ρ_DM.

A more meaningful comparison can be made between the cumulative total mass profiles of the two clusters as obtained from our merger simulation XDBf_sb, against those found by the lensing reconstruction of the El Gordo cluster mass profiles. For the NW and SE cluster Fig. 18 of Kim et al. (2021) shows separately the total halo mass with a given radius, these profiles can be contrasted with the corresponding ones extracted from simulation XDBf_sb and plotted in Fig. 14. A visual comparison between the two figures overall shows that for each cluster the corresponding mass profiles exhibit approximately the same radial behavior. However, at large radii the mass profile of the NW cluster recovered from lensing data is not described well by the corresponding profile extracted from the simulation, which tends to be systematically higher by about a factor of 1.5 − 2.

To be more quantitative in each panel we quote the values of M₂₀₀ and r₂₀₀, as determined according to Eq. (3) from the cumulative mass profile of the cluster at the redshift z = 0.87. These values can be compared with the corresponding mass estimates of Kim et al. (2021): $M_{200}^{\mathrm{NW}}(K)=9.9_{-2.2}^{+2.1}\times {10}^{14}{M}_{\odot }$ and $M_{200}^{\mathrm{SE}}(K)=6.5_{-1.4}^{+1.9}\times {10}^{14}{M}_{\odot }$ for the NW and SE cluster, respectively.

From the right panel of Fig. 14 it can be seen that the mass $M_{200}^{\mathrm{NW}}(t=0.24\mathrm{Gyr})\sim 2\times {10}^{15}{M}_{\odot }$ of the NW cluster is a factor of ∼2 higher than the corresponding lensing mass estimate $M_{200}^{\mathrm{NW}}(K)$. This is not surprising, as the initial mass of the primary of model XDB_sb was set to $M_{200}^{(1)}=1.6\times {10}^{15}{M}_{\odot }$ because its collision parameters are the same of model Bf_rc20. According to the findings of Sect. 3.1.2, this is the fiducial models that best matches the observations. The collision parameters of this model (see Table 7) have been therefore specifically chosen to perform the SIDM merger simulations of Fig. 12.

For the SE cluster we find $M_{200}^{\mathrm{SE}}(t=0.24\mathrm{Gyr})\sim 9\times {10}^{14}{M}_{\odot }$, marginally inconsistent at 1σ level with the estimate $M_{200}^{\mathrm{SE}}(K)$. This mild tension is easily understood to originate from the assumed initial masses M₂₀₀ for the two colliding clusters. From Fig. 14 the present cluster masses at r = r₂₀₀ are found to be higher by a factor of ∼30% compared to their initial values. This is due to the flattening of the DM inner density profiles during the merger, and at the present epoch this in turn leads to an average higher DM density at large cluster radii. For the SE cluster initially $M_{200}^{\mathrm{SE}}=6.5\times {10}^{14}{M}_{\odot }$ and r₂₀₀ ∼ 1.3 Mpc; however due the presence of cored DM profiles after the core passage, Eq. (3) is now solved for $M_{200}^{\mathrm{SE}}\sim 8.9\times {10}^{14}{M}_{\odot }$ and r₂₀₀ ∼ 1.43 Mpc.

These discrepancies between the mass profiles tend to become less severe at inner radii. We chose R = 500 kpc as the minimum cluster radius at which the cluster mass profiles of Fig. 18 (Kim et al. 2021) start to appreciably diverge. The masses of the NW and SE cluster at this radius are then estimated to be about $M_{\mathrm{extr}}^{\mathrm{NW}}(R<500\text{kpc})\simeq 3\times {10}^{14}{M}_{\odot }$ and $M_{\mathrm{extr}}^{\mathrm{SE}}(R<500\text{kpc})\simeq 2.34\times {10}^{14}{M}_{\odot }$, respectively.

These masses are indicated in the corresponding panels of Fig. 14 with magenta triangles. For the NW cluster the difference between the cluster mass $M_{\mathrm{sim}}^{\mathrm{NW}}(R<500\text{kpc})$, as extrapolated from the mass profile of the simulation, and $M_{\mathrm{extr}}^{\mathrm{NW}}$ is still significant. However, for the SE cluster the accord between $M_{\mathrm{sim}}^{\mathrm{SE}}(R<500\text{kpc})$ and $M_{\mathrm{extr}}^{\mathrm{SE}}$ is now much better.

We then conclude that present data do not allow us to firmly rule out the presence of cored DM profiles in the El Gordo cluster. Moreover, we argue that these differences between the predicted and reconstructed mass profiles can be interpreted in terms of the adopted initial conditions for model XDB_sb. For instance, we suggest that setting the initial mass of the secondary for the merger model XDB_sb to approximately ∼20% less than the initial adopted value of $M_{200}^{(2)}\sim 6.5\times {10}^{14}{M}_{\odot }$ will provide, for the SE cluster, a mass profile in better agreement with observational estimates.

However, it must be stressed that WL-based cluster mass measurements (see Umetsu 2020, and references cited therein) are subject to a number of uncertainties. To break the mass-sheet degeneracy (Umetsu 2020) cluster masses must be estimated by assuming a specific mass profile, such as an NFW. In a recent paper, Lee et al. (2023) showed that in major cluster mergers (∼10¹⁵ M_⊙) departure from the assumed mass profile can be significant and the subsequent WL mass bias can be as high as ∼50%. They conclude that in massive mergers previous WL mass estimates can be significantly overestimated.

Nonetheless, we suggest that in the case of El Gordo cluster such a large range of mass uncertainties should be taken with some caution. This is motivated by recent, independently based, SL mass estimates. Kim et al. (2021) found a WL total mass estimate of ∼2 × 10¹⁵ M_⊙ for the El Gordo cluster, this is roughly in the same range as the corresponding mass estimates derived for El Gordo cluster from recent SL analyses (Caminha et al. 2023; Diego et al. 2023).

3.4. Testing the dependence of the X-ray to DM peak offset on the initial halo concentration and gas density profile of the SE cluster

As we have seen in Sect. 3.2, a returning scenario fails in reproducing the observed displacement $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ of the X-ray peak from the SE mass centroid. On the contrary, the results of the previous Section demonstrate that a collisional DM scenario can naturally account for the observed offset. The size at the observer epoch of the observed offset $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ depends however on a number of additional factors that can complicate the interpretation of $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ in terms of an SIDM scenario.

The first physical effect during the merger that can have a significant impact on the motion of the cool core of the secondary is the depth of the potential well of the hosting DM halo. An initially shallower inner DM density profile will in turn lead to a less deep potential in the cluster inner regions, thereby easing the motion of the SE cool-core during the cluster merger.

In order to assess the impact of a shallower DM potential on the final displacement of the SE cluster cool core, in this section we present results from merger simulations with initial condition setups identical to those of previous runs, but with different values of the initial halo concentration parameter ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}$ of the SE cluster. To this end, we use as reference the merger models DBf_rc20 of Fig. 10 and XDBf_sb of Fig. 12. For both of these models, for the mass of the secondary one has from Tables 6 and 7 $M_{200}^{\mathrm{SE}}=6.9\times {10}^{14}{M}_{\odot }$. According to Eq. (5), for these simulations the initial halo concentration parameter is then set to ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2.682$.

For each of these two merger models, we ran two alternative simulations. The initial condition setup of these simulations is identical to that of the parent simulation; the only difference is in the initial value of the halo concentration parameter for the SE cluster. For the two simulations we set this parameter to ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2\equiv c_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and to ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=3.2\equiv c_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}$, respectively. This interval of values for the concentration parameter approximately corresponds to what is predicted by the standard deviation (σ_c/c ∼ 0.3) of the M − c relation (Bhattacharya et al. 2014), and it was specifically chosen in order to bracket, for different shapes of the DM potential, the motion of the SE cool core during the merger.

Figure 15 shows the X-ray maps, as extracted at the observer epoch from our alternative merger simulations. We first discuss results from the two alternative simulations (top panels) of the parent merger model DBf_rc20. A striking feature that emerges from the simulation results is that the final size of the $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ offset is heavily influenced by the shape of the density profile of the hosting DM halo. The top-left panel shows that for the DBf_rc20 run with ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2$ the X-ray peak is trailing the DM centroid, with the relative offset now negative ($d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}\simeq -190\text{kpc}$).

This offset between the X-ray peak and the SE mass centroid was already present in the merger models DBf_rc20 of Fig. 10, as expected because of the different behavior during the merger between the dissipative ICM and the collisionless DM component. However, the size of $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ is now much larger, thus showing that for the SE cluster the motion of the cool core during the merger is very sensitive to the choice of relatively low values of c₂₀₀. On the contrary, the top right panel of Fig. 15 shows that when ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=3.2$ the size of $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ is unaffected by the choice of c₂₀₀ and is very similar to that of the standard run.

Model XDBf_sb is an SIDM merger model with σ_DM/m_X = 4 cm² gr⁻¹, the bottom-left and bottom-right panels of Fig. 15 show the X-ray maps of the corresponding ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2$ and ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=3.2$ simulations, respectively. An important feature that is clearly seen from the bottom-left panel of Fig. 15 is that for the ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2$ simulation of model XDBf_sb the $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ offset is now much smaller than in the corresponding DBf_rc20 run.

This behavior is the consequence of two competing effects: the easier motion of the cool core during the merger because of the shallower DM density profile (as in the BDf_rc20), and the deceleration of the DM halos (as in XDBf_sb) due to the drag force induced during the cluster collision by the scatterings of the DM particles, with the subsequent loss of momentum.

From the bottom-right panel of Fig. 15, it can be seen that for the ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=3.2$ simulation of model XDBf_sb the $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ offset is insensitive to the choice of c₂₀₀, as in the corresponding DBf_rc20 run. We argue that for values of ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}\gtrsim 2.7=c_{d_{\mathrm{X}-\mathrm{D}\mathrm{M}}}$, the size of the final offset between the X-ray and DM peak is weakly sensitive to the form of the inner DM potential of the SE cluster.

We also conclude that these findings are not affected by observational constraints on the range of allowed values for ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}$. As previously discussed, from their WL study on El Gordo cluster Kim et al. (2021) obtained ${\text{c}}_{200}^{\text{SE}}(K)=3.2_{-0.74}^{+1.17}$ for the posteriors of the concentration parameter of the SE cluster. Their Fig. 7 (bottom panel) shows that their estimate ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}(K)\sim 3.2$ coincides approximately with c_high. This shows that the 1σ lower limit on ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}(K)$ and the previously derived threshold c_{dX − DM} have approximately the same value, thus indicating that the size of the offset $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ is not affected by observational uncertainties for ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}$.

The size of the X-ray cool core is another important factor that contributes to its final displacement from the SE mass peak. As outlined in the Introduction, the simplest physical interpretation of the El Gordo merger (Menanteau et al. 2012; Jee et al. 2014) is a scenario in which the original gas core of the primary has been destroyed by the collision with the compact core of the secondary, during its motion from NW to SE. Clearly, a very dense, compact, cool core will be less prone to the ram pressure experienced during its motion through the ICM of the primary.

Following ZYL15, for the initial gas density profile of the secondary we adopted the Burkert (1995) density profile (Eq. (9)), and set the gas core radius to r_c = r_s/3. This is a reasonable choice that has been found to reproduced fairly well the observed X-ray luminosity L_X (see Table 4), but there are other choices as well.

We want then to study the impact on the final position of the SE cluster X-ray peak of initially different gas density profiles. To this end we used the merger model XDBf_sb as reference and we performed simulations with the same initial condition setup but with different initial gas density profiles for the SE cluster. In particular, we chose to model the gas density profile of the secondary according to the modified beta profile of Vikhlinin et al. (2006):

$$\begin{array}{c}\hfill {\rho }_g^2\propto \frac{{(r/\underset{c}{\overset{\sim }{r0pt0.8ex}})}^{-\alpha }}{{[1+{(r/\underset{c}{\overset{\sim }{r0pt0.8ex}})}^2]}^{3\beta -\alpha /2}}\frac{1}{{[1+{(r/\underset{s}{\overset{\sim }{r0pt0.8ex}})}^{\gamma }]}^{ϵ/\gamma }}·\end{array}$$(28)

The shape of this profiles is controlled by six parameters: $\{\alpha ,\beta ,\gamma ,ϵ,{\tilde{r}}_c,{\tilde{r}}_s\}$, thus allowing a large amount of freedom in regulating its radial behavior. Following Chadayammuri et al. (2022), we set here β = 2/3, γ = 3 and ϵ = 3. The size of the cool core is then controlled by ${\tilde{r}}_c$ and α, while the scale radius ${\tilde{r}}_s$ regulates the transition to ∝r^{−(3β + ϵ/2)}.

To ease the comparison with previous findings (Chadayammuri et al. 2022), it is convenient to express the radii ${\tilde{r}}_c$ and ${\tilde{r}}_s$ in units of a_s, the scale radius of the generalized “super-NFW” (sNFW) profile (Lilley et al. 2018):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}(r)=\frac{3M}{16\pi a_s^3}\frac{1}{(r/a_s){(1+r/a_s)}^{5/2}}·\end{array}$$(29)

The scale radius, a_s, can be recovered from the NFW concentration parameter c_NFW according to the linear relation c_sNFW = 1.36 + 0.76c_NFW, where c_sNFW ≡ 3r₂₀₀/2a_s is the sNFW concentration parameter. For the SE cluster of model XDBf_sb one has $r_{200}^{\mathrm{SE}}\sim 1.32\text{Mpc}$, ${\mathrm{c}}_{200}^{\mathrm{S}\mathrm{E}}=2.68$ and a_s ∼ 0.58 Mpc.

Figure 16 presents for different values of the profile parameters $\{{\tilde{r}}_c,{\tilde{r}}_s,\alpha \}$ a set of initial gas density profiles, as described by Eq. (28). For the sake of comparison we also show the Burkert profile (Eq. (9)), previously adopted for the SE cluster in all of the merger models. In the panel, the radii ${\tilde{r}}_c$ and ${\tilde{r}}_s$ are expressed in units of a_s and all of the profiles are consistently normalized to $f_g^{\mathrm{SE}}=0.14$ at r = r₂₀₀.

As can be seen, the size of inner cool core region and the cuspiness of the central gas density peak can be regulated by an appropriate choice of the profile parameters and range from that of a compact cool core (model Vk_a: $\{{\tilde{r}}_c=0.05,{\tilde{r}}_s=0.6,\alpha =2\}$) to that of an almost non-cool core (model Vk_e: $\{{\tilde{r}}_c=0.4,{\tilde{r}}_s=0.9,\alpha =0.5\}$). This dependence of the radial profile behavior on the set of parameters $\{{\tilde{r}}_c,{\tilde{r}}_s,\alpha \}$ can be compared with similar plots shown in Fig. 12 (left panel) of Chadayammuri et al. (2022).

To explore the dependence of the observed offset $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ on the initial gas density profile of the SE cluster, as for the merger models of Fig. 15 we ran two alternative simulations of model XDBf_sb. The initial condition setup of these simulations is identical to that of the parent simulation, the only difference being now in the choice of initial gas density profile for the SE cluster.

In order to bracket a plausible range for the compactness of the cool core, for the two alternative simulations we chose to model the gas density profile of the secondary according to models Vk_b and Vk_e of Fig. 16, respectively. The X-ray maps extracted from these two alternative merger simulations are shown in Fig. 17.

For model Vk_b the left panel of Fig. 17 shows a $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ offset nearly identical ($d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}\sim 90\text{kpc}$) to that of the parent simulation XDBf_sb of Fig. 12. This demonstrates that, above a certain threshold, the size of the final offset of the X-ray peak with respect to the DM centroid can be considered insensitive to the cuspiness of the SE cluster inner gas density. On the other hand, it worth noting that for this model the X-ray luminosity (L_X ≃ 7 × 10⁴⁵ ergs⁻¹) is now higher by a factor of ∼3 with respect the corresponding observational range as reported in Table 4.

For model Vk_e (right panel of Fig. 17), the $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ offset is now negative ($d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}\sim -60\text{kpc}$), with the X-ray peak now trailing the DM centroid. We interpret this behavior as a consequence of the much less compact gas density peak with respect that of model Vk_b, this implying a larger ram pressure force, and in turn a larger deceleration, experienced by the secondary’s cool core during the cluster collision (Ricker & Sarazin 2001; Poole et al. 2006, 2008; Mitchell et al. 2009).

To summarize, our results demonstrate that both the halo concentration parameter and the initial shape of the gas density profile of the SE cluster have an impact on the offset $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ of the X-ray peak relative to that of the DM. We found that both of these dependences could act to reduce the offset size, but are however unable to increase it above the observational estimate.

Our conclusion is therefore that the interpretation of the observed offset between the X-ray peak and the DM centroid in an SIDM scenario remains valid. The alternative configurations of the initial condition setup analyzed here for the SE cluster are found to be ineffective in increasing the size of the offset.

Finally, it can be seen from Fig. 17 that very compact cool core gas density profiles can produce $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$ offsets in good accord with the measured value, but they fail to account for the observed X-ray luminosity. This strongly suggests that the previous used Burkert profile (Eq. (9)), with the gas core radius set to r_c = r_s/3, can be considered a good approximation to the gas density profile of the SE cluster.

4. Summary and conclusions

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

To this end, we considered a wide range of initial conditions, performing both off-axis and head-on merging runs. The initial conditions of the binary merging clusters were constrained to be consistent with recent lensing-based mass measurements (Diego et al. 2020; Kim et al. 2021). We imposed the following range of values for the mass of the primary: 10¹⁵ M_⊙ ≲ M ≲ 1.6 × 10¹⁵ M_⊙. A summary of our main results is as follows:

(i) The observed twin-tailed X-ray morphology, as well as other observational constraints, are well matched by the off-center fiducial merger models Bf_rc20, Bg_rc20, and Bl_rc24 (Sect. 3.1), which have collision velocities and impact parameters in the range 2000 km s⁻¹ ≲ V ≲ 2500 km s⁻¹ and 600 kpc ≲ P ≲ 800 kpc, respectively.

These findings demonstrate that in simulations with cluster masses lower than previously considered (Donnert 2014; Molnar & Broadhurst 2015; Zhang et al. 2015, 2018), it is possible to consistently reproduce the observed twin-tailed X-ray morphology seen in the cluster El Gordo. Finally, the observational constraints discussed in Sect. 3.1 tend to favor merger models with a primary mass of about M ≃ 1.6 × 10¹⁵ M_⊙.

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

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

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

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

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

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

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

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

The most important striking feature that emerges from the SIDM merging models of the El Gordo cluster is that the observed offsets between the different mass components are well reproduced by the simulations XDBf_sa and XDBf_sb (Table 7). These simulations have a DM cross-section in the range σ_DM/m_X ∼ 4 − 5 cm² gr⁻¹ and initial collision parameters of model Bf_rc20 (Table 2).

Another significant aspect of these merger models is the value of the mean relative line-of-sight radial velocity between the two BCGs, which in contrast with previous runs (see Fig. 7) is now on the order of several hundred km s⁻¹ and in much better agreement with data. As discussed in Sect. 3.3, this is a direct consequence of the gravitational pull experienced by the two BCGs because of the dissipative interactions now present between the two DM halos.

Although these models are not able to reproduce the observed X-ray morphology as well as in the standard CDM runs (Sect. 3.1), we nonetheless argue that these findings provide significant support for the possibility that the DM behavior in the El Gordo cluster exhibits collisional properties. Additionally, it must be stressed that a significant advantage of the SIDM merging models is that they offer a natural explanation for all of the observed offsets, while in the standard CDM scenario the observed separations between the different centroids remain either unexplained (X-ray peak to DM) or fraught with difficulties (BCG to DM).

From these conclusions follow two main consequences that need to be discussed. The first is that the two SIDM models of Table 7 do not admit free parameters; the observed offsets are well reproduced only for a DM scattering cross-section in the range σ_DM/m_X ∼ 4 − 5 cm² gr⁻¹ and initial collision parameters given by model Bf_rc20. This is a direct result of the post-pericenter evolution in an SIDM scenario, for which the observed separations between the different centroids are locked to the measured distance, d_DM, between the DM mass centroids, which must be about d_DM ∼ 700 kpc at the present epoch.

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

The other striking feature of the SIDM merger models XDBf_sa and XDBf_sb is that the best match to the data is obtained only for values of σ_DM/m_X in the very narrow range between ∼4 cm² gr⁻¹ and ∼5 cm² gr⁻¹. Such values are strongly excluded by the present upper limits on the SIDM cross-section, which on galaxy cluster scales have been derived independently by various authors (Robertson et al. 2017a, 2019; Kim et al. 2017; Wittman et al. 2018; Harvey et al. 2019; Shen et al. 2022; Cross et al. 2024). These upper limits are at best of order unity (σ_DM/m_X ≲ 1 cm² gr⁻¹), although Wittman et al. (2018) claim that the constraint can be relaxed up to σ_DM/m_X ≲ 2 cm² gr⁻¹ by properly averaging galaxy-DM offsets in a sample of merging clusters.

For example, in the case of the Bullet Cluster, the absence, within the observational errors on the centroid positions, of an observed galaxy-DM separation (≲20 kpc, Randall et al. 2008) has been used in SIDM merging simulations (Markevitch et al. 2004; Randall et al. 2008; Robertson et al. 2017a) to constrain the allowed values of the DM cross-section to lie below σ_DM/m_X ≲ 1 cm² gr⁻¹. This evidently implies that an SIDM merging simulation of the Bullet Cluster, with the values derived here for the DM cross-section (σ_DM/m_X ≲ 4 − 5 cm² gr⁻¹), would lead to simulation results clearly incompatible with the observational constraints.

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

Kim et al. (2017) present a set of SIDM merging simulations of equal-mass massive clusters, with masses on the order of ∼10¹⁵ M_⊙. They compared the expected galaxy-DM offsets from the SIDM simulations against measured galaxy-DM offsets in observed mergers. The authors were able to derive an upper limit on the galaxy-DM offset of ≲20 kpc for σ_DM/m_X = 1 cm² gr⁻¹, in accordance with the range of the spatial separations shown in the left panels of Fig. 11 for the σ_DM/m_X = 1 cm² gr⁻¹ simulations. As noticed by the authors, these bounds are an order of magnitude smaller than the measured offsets for these merging systems (∼100 − 300 kpc), and thus in tension with their simulations. Their conclusions were, however, affected by the large observational uncertainties for the measured offsets in massive mergers of ∼10¹⁵ M_⊙.

Nonetheless, according to Kim et al. (2021), the presence of offsets between the different mass components of the El Gordo cluster is statistically significant, which in turn implies that the null hypothesis of zero size offsets is statistically disfavored. On the other hand, is difficult to assess the statistical significance of the offset sizes we found for the merger model XDBf_sb (Fig. 12) with σ_DM/m_X = 4 cm² gr⁻¹. This indeterminacy is caused by the absence of positional errors for the various peak locations used to calculate the offsets.

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

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

The size of positional errors for the X-ray emission peak of the SE cluster is expected to be less significant than for the other centroids, due to the squared dependence of the X-ray emission with the gas density. These errors are quantified according to the observational procedure adopted to process the raw X-ray images and to extract the positions of substructures in the X-ray emission, together with the relative errors. In the case of El Gordo, the effective positional error of the X-ray peak is relatively small and set by the angular resolution (∼0.5″) of the Chandra X-ray image².

The uncertainty in the size of the observed offset between the X-ray and the SE mass peak is then dominated by the uncertainty σ_DM in the position of the mass peak; this allows us to place the offset in the range $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}\sim 100\pm 40\text{kpc}$. An important result that follows from this finding is that SIDM merger simulations of the El Gordo cluster with σ_DM/m_X ≲ 2 cm² gr⁻¹ are marginally inconsistent with the observed offset, $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}$. As can be seen from the maps of Fig. 11, for these models the offset of the X-ray peak is always limited to $d_{\mathrm{X}-\mathrm{DM}}^{\mathrm{SE}}\lesssim 40\text{kpc}$. In the following we thus adopt the view that the estimated offset sizes for the El Gordo cluster are in the closest agreement with SIDM merger simulations for values of the DM cross-section in the range σ_DM/m_X ∼ 4 − 5 cm² gr⁻¹.

To resolve the inconsistency between this range of values for σ_DM/m_X and its current upper limits on cluster scales (σ_DM/m_X ≲ 1 cm² gr⁻¹), we suggest that the SIDM model proposed so far should only be considered a first-order approximation to the modeling of the physical processes describing DM interaction. Specifically, we argue that the capacity of DM to exhibit collisional properties during a cluster merger is correlated with some other feature of the collision itself.

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

It must be stressed that such an energy behavior is completely different from the mechanism previously proposed (Kaplinghat et al. 2016), in which a power-law dependence of σ_DM/m_X on the collision velocity has been postulated to satisfy its present upper limits over a wide range of scales. Such a consideration also applies to other SIDM models (Tulin & Yu 2018) in which DM particles scatter inelastically (Vogelsberger et al. 2019) or anisotropically (Robertson et al. 2017b).

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

To validate the proposed picture, we took as reference the El Gordo cluster, for which the measured galaxy-DM and BCG-DM offsets lie in the range 50 − 100 kpc. The best match to the data is obtained for models XDBf_sa and XDBf_sb, which have the initial collision parameters $\{M_{200}^{(1)},q,V,d_{\mathrm{ini}}\}$ of model Bf_rc20: {1.6 × 10¹⁵ M_⊙,  2.3,  2500 km s⁻¹,  6.1 Mpc}, where $d_{\mathrm{ini}}=2(r_{200}^1+r_{200}^2)\simeq 2(1.74+1.32)$. With these parameters we can estimate using Eq. (17) the energy of the collision: E_EG ∼ 1.4 × 10⁶⁴ ergs.

We then compared this finding to a merging system for which no appreciable offsets have been detected between the DM and the stellar component: the Bullet Cluster. This is a well-known cluster merger that has been the subject of many numerical studies (Springel & Farrar 2007; Mastropietro & Burkert 2008; Lage & Farrar 2014, 2015; Robertson et al. 2017a). We extracted the collision parameters $\{M_{200}^{(1)},q,V,d_{\mathrm{ini}}\}$ from Table 1 (LF entry) of Lage & Farrar (2015): {1.9 × 10¹⁵ M_⊙,  7.3,  2800 km s⁻¹,  2.8 Mpc}. According to these collision parameters, the total energy content (17) of the Bullet Cluster is E_Bullet ∼ 3 × 10⁶³ ergs.

According to our hypothesis, these simple estimates bracket the allowed values of E_crit to be approximately in the range E_Bullet ≲ E_crit ≲ E_EG. To further corroborate these findings, it would be interesting to verify if there are other major merger clusters that exhibit large galaxy-DM peak offsets (∼100 kpc). If the merging clusters are as massive as El Gordo, this would clearly lend support to the idea that there is a critical energy threshold, E_crit, around the value of E_EG.

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

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

We estimate the collision parameters of the Sausage Cluster from Jee et al. (2015): $\{M_{200}^{(1)},q,V,d_{\mathrm{ini}}\}=\{1.1\times {10}^{15}{M}_{\odot },1.12,1300\text{km}{\text{s}}^{-1},4\text{Mpc}\}$, where $d_{\mathrm{ini}}=2(r_{200}^1+r_{200}^2)\sim 4\text{Mpc}$ and V is estimated from Eq. (14) of Jee et al. (2015) by setting b = 0 and t_impact ∼ 11 Gyr. We can now find the energy of the merging, which turns out to be about E_Sausage ∼ 1.5 × 10⁶⁴ ergs.

This energy value is very close to that of E_EG and strongly supports the hypothesis of the existence of an energy threshold that regulates DM behavior. This implies that, during the cluster collision, DM can exhibit collisional properties according to the amount of collisional energy.

Numerical simulations of the Sausage Cluster have previously been performed (van Weeren et al. 2011; Donnert et al. 2017; Molnar & Broadhurst 2017), but without including the effect of SIDM. From the SIDM merging simulations of Kim et al. (2017), the expected galaxy-DM offset for merging clusters of ∼10¹⁵ M_⊙ is on the order of ≲20 kpc for σ_DM/m_X = 1 cm² gr⁻¹, which is smaller than the estimated offset for the Sausage Cluster by about a factor of 10.

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

---

Acknowledgments

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

References

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

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

The Monte Carlo method is then implemented as follows. We assumed that the parallel code subdivides the computation into P tasks, and the simulation volume is correspondingly split into Ω_p non-overlapping sub-volumes, with p = 0, 1, …, P − 1 being the task index. The load balancing domain decomposition is done by cutting segments of a globally ordered Hilbert space-filling curve, as in the code GADGET-2 (Springel 2005).

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

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

To avoid multiple scatterings, it is now necessary that each task p knows if a member of a scattering pair $\{i,j{\}}_p^{oth}$ has undergone a scattering event in another task, this in order to accordingly update its velocity. To accomplish this for each pair $\{i,j{\}}_p^{oth}$ we associate a 64 bit integer key defined as $I{D}_p^{(i_K,j_K)}=i_K(i)+M{j}_K(j)$, where i_K(i) denotes an integer key of the local DM particle i. The integer key i_K(i) is unique for each DM particle and it is assigned initially; it is a DM particle property that moves with the particle when this migrates between processors. The integer M (2³¹) must be significantly larger than N_DM; therefore, the integers $I{D}_p^{(i_K,j_K)}$ are unique for each scattering pair in the P tasks.

To consistently incorporate multiple scatterings in the parallel code, we now perform a global sorting of the P sets of integers $I{D}_p^{(i_K,j_K)}$, with p = 0, 1, …, P − 1. After the sorting each task p knows the global list of sorted keys, and therefore all of the other scatterings and subsequently in which task they are taking place.

Finally, the post-scattering velocities (20) of each scattering pair in task p are evaluated as follows. After the global sorting each task will process its local set of scattering pairs by looping iteratively over its local list of sorted keys, this in increasing key order because of the sorting. Pairs that exhibit a scattering dependence with $I{D}_p^{(i_K,j_K)}>I{D}_q^{(j_K,m_K)}$ are put in a waiting list until the updated velocities of particle j have been received from task q. We note that this implementation allows for a single particle to undergo multiple scatterings.

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

To test the SIDM part of the code we evolve an isolated DM NFW halo, initially at equilibrium, for different values of the SIDM cross-section. We considered the following values: σ_DM/m_X = 0,  1, and 10cm² gr⁻¹. For the DM halo, we set M₂₀₀ = 1.3 ⋅ 10¹⁵M_⊙ at z = 0.87, so R₂₀₀ = 1.63 Mpc and for the concentration parameter c₂₀₀ = 2.55 from Equation (5). The initial condition setup of the DM halo is described in Section 2.2.1.

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

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathit{th}}(r)=\frac{{\sigma }_{\mathit{DM}}}{m_X}<v_{\mathit{pair}}(r)>{\rho }_{\mathit{DM}}(r),\end{array}$$(A.1)

where ρ_DM(r) is the local DM density and < v_pair(r)> the mean pairwise velocity. Assuming a Maxwell-Boltzmann distribution this quantity can be written as (Robertson et al. 2015) $<v_{\mathit{pair}}(r)>=(4/\sqrt{\pi }){\sigma }_{1D}(r)$, where σ_1D(r) is the 1D velocity dispersion of the NFW halo (Łokas & Mamon 2001).

To contrast the predicted rate, Γ_th(r), with the corresponding rate extracted from the simulations, Γ_sim(r), we measured Γ_sim(r) as follows. We first accumulated the number of scattering events as a function of radius in radial bins together with the number of particles in the same bin. At epoch t, the binned scattering rate, Γ_sim(r), was then obtained by dividing for each bin the number of scattering events by the number of particles. It must be stressed that during the simulations the rates Γ_sim(r) are obtained by evaluating the local scattering probability (19), but without updating the particle velocities (20) of the pair. This is to avoid any modification to the DM density profile, which otherwise will make Equation (A.1) not applicable.

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

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

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

Following Springel et al. (2005) we relate the scale radius r_s and concentration parameter c of the NFW profile to the scale length a_H of an Hernquist density profile with the same mass within r₂₀₀:

$$\begin{array}{c}\hfill a_H=r_s\sqrt{2[ln(1+c)-c/(1+c)]}.\end{array}$$(A.2)

For the halo parameters we used r_s = r₂₀₀/c ≃ 0.64Mpc and a_H ≃ 660kpc, which is only ∼10% higher than $a_H^Z$. Thus, at inner radii our profiles can be compared with the DM halo profiles displayed in top-left panel of Figure 1 of ZuHone et al. (2019). From the left panel of Figure A.1 at the initial epoch (t = 0) we estimate at the radius r ∼ 10⁻²r₂₀₀ ≃ 16kpc a central DM density of ρ_DM/ρ_c ≃ 8 ⋅ 10⁴. For the chosen redshift (z = 0.87) this value corresponds to ρ_DM ≃ 3 ⋅ 10⁷M_⊙kpc⁻³, in accordance with the corresponding value estimated from Figure 1 of ZuHone et al. (2019).

We now estimate at t = 1 Gyr the central value of the halo DM density for the SIDM run with σ_DM/m_X = 10cm² gr⁻¹. For this test from Figure A.1 we obtain at the radius r ∼ 10⁻²r₂₀₀ a central DM density of ρ_DM/ρ_c ≃ 1.8 ⋅ 10⁴, which is equivalent to ρ_DM ≃ 6.6 ⋅ 10⁶M_⊙kpc⁻³. For the same SIDM run this value is in close agreement with the corresponding DM density, estimated at the same radius and at the same epoch from Figure 1 of ZuHone et al. (2019). The results of our tests then validate our implementation of DM self-interactions, and demonstrate that the resulting code can be safely used to perform the SIDM simulations presented here.

All Tables

All Figures