© The Authors 2025

1. Introduction

According to the standard hierarchical structure formation model, larger cosmic structures are formed through the merging of smaller structures. Galaxy clusters are created at the intersections of the cosmic web through the merging process of smaller substructures. In the later stages of structure formation, nearby galaxy clusters start merging with each other because of their mutual gravitational attraction. Galaxy cluster mergers are the most energetic phenomena in the Universe, with colliding velocities of ∼1000 km s⁻¹ and a released energy of 10⁶⁴ ergs (Sarazin 2002). Cluster mergers offer a unique opportunity to study numerous physical processes ranging from the acceleration of particles and cosmic rays to cosmology. For instance, valuable insights on the self-interaction of dark matter particles can be obtained by studying cluster mergers, particularly the offset between the dark matter and gas distributions in post-mergers (Markevitch et al. 2004; Robertson et al. 2017). Furthermore, the kinetic energy of the clusters can dissipate to the intracluster medium (ICM) and produce shock and cold fronts detected in X-rays by their density and temperature jumps. Such shock and cold fronts can provide information on the microphysics of the ICM (e.g., thermal conduction and magnetic fields, Markevitch et al. 2003; Chadayammuri et al. 2022; Sarkar et al. 2023) and its energy transport processes as well as on the dynamics of merging clusters.

On the other hand, galaxy cluster pre-mergers constitute ideal targets for studying gas bridges that connect individual halos. Analyzing such gas bridges can help provide better understanding of the location of the missing baryons in the late Universe (e.g., Werner et al. 2008). Shock fronts were also shown to exist in gas bridges between pre-merging clusters (e.g., Kato et al. 2015; Sarkar et al. 2022; Omiya et al. 2023) due to colliding gas. Because of the great interest in them, galaxy cluster mergers have been extensively studied in the last decades in the full electromagnetic spectrum and at all merging stages (e.g., Markevitch et al. 2005; Bradač et al. 2006; Cassano et al. 2010; van Weeren et al. 2010; Tempel et al. 2017; Ha et al. 2018).

Moreover, the study of the individual clusters comprising pre-merging systems can help determine if they can be used as typical clusters for cosmological and astrophysical studies, such as for constraining the halo mass function and galaxy cluster scaling relations. Interacting clusters are usually disturbed systems, with their level of relaxation (as determined from, e.g., their morphological parameters) depending on their merging stage. Disturbed clusters have been shown to often follow scaling relations different from those of the general cluster population (e.g., Lovisari et al. 2020; Migkas et al. 2021), which can bias their use as cosmological probes. However, the studies showing this include the entire systems in their analysis (i.e., merging and non-merging cluster sides). No extensive effort has been made to characterize only the non-interacting sides of pre-merging systems and determine if their behavior is consistent with isolated relaxed clusters.

X-ray observations are especially useful for studying cluster mergers. X-ray data from instruments such as XMM-Newton, Chandra, and eROSITA can spatially resolve the hot gas of merging clusters more efficiently than Sunyaev-Zeldovich (SZ) observations due to their superior angular resolution. X-ray observations also allow for the study of the properties of the ICM in much greater detail. Additionally, optical data trace the galaxy population of clusters and their weak and strong lensing signal. Combining X-ray and optical information can provide a complete picture of the state of a merging cluster system and its components.

In this work, we perform the first X-ray analysis of a previously identified cluster merger, SPT-CLJ2228-5828, alongside a refined weak lensing (WL) analysis of its cluster components. Based on earlier data, the cluster was a strong candidate for a post-merger system (see Sect. 2). However, we find that SPT-CLJ2228-5828 consists of two individual clusters in a pre-merging phase. The collision of the gas of the two clusters has started, with a gas bridge between the components where a shock front can be found. While the two clusters are highly similar in X-rays, only one of these clusters shows a clear WL detection, with the other being marginally detected.

The structure of the paper is as follows: In Sect. 2, we present the results of past analyses of the SPT-CLJ2228-5828 merging system. In Sect. 3, we describe the X-ray and optical data used in this work, the data reduction, the imaging, the X-ray spectral analysis, and the WL analysis. In Sect. 4, we present the results regarding the two cluster members of the merger. In Sect. 5, we present the analysis of the gas bridge connecting the two clusters. Finally, in Sect. 6, we discuss and summarize our findings. Throughout the paper, a ΛCDM model is assumed with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3.

2. The SPT-CLJ2228-5828 cluster system

The SPT-CLJ2228-5828 system was first detected as a single cluster in the 2500 deg² SPT-SZ survey by Bleem et al. (2015) at (RA, Dec) = (337.215° , − 58.469°) and a photometric redshift z = 0.71 ± 0.06. From its integrated total Compton parameter the system’s total mass within R₅₀₀ was estimated to be M₅₀₀ = (3.17 ± 0.82)×10¹⁴ M_⊙¹, pointing to a cluster of average mass. Its apparent radius was determined to be 1.5′, corresponding to ≈650 kpc and classifying it as one of the 43 most disturbed clusters from the 2500 deg² SPT-SZ survey sample (Bleem et al. 2015). A follow-up analysis of a subsample of the 2500 deg² SPT-SZ catalog by Bocquet et al. (2019) yielded similar results with $M_{500}=3.{27}_{-0.83}^{+0.63}\times {10}^{14}$ M_⊙ and z = 0.734 ± 0.046 (photo-z), with the system still considered a single cluster. SPT-CLJ2228-5828 was also detected as a single object by Hilton et al. (2021) using the Atacama Cosmology Telescope and the SZ effect. The determined redshift and mass of the system were z = 0.760 ± 0.013 (photo-z) and M₅₀₀ = (3.82 ± 0.79)×10¹⁴ respectively, consistent with previous studies. Zenteno et al. (2020) performed a joint SZ and optical analysis of the system that revealed two distinct galaxy overdensities, with the cluster gas, as traced by the SPT-SZ data, located in the middle of these two. The identified brightest cluster galaxy (BCG) of the two galaxy distributions was found to be separated by 0.49 R₂₀₀ from the SZ centroid, corresponding to 1.4′ and ≈605 kpc. Furthermore, Schrabback et al. (2021) performed a WL analysis of the system using shallow Hubble Space Telescope (HST) data and found that there are two peaks in the WL mass reconstruction ≈1 Mpc apart from each other, coinciding with the two galaxy overdensities. Based on the results of the two referenced studies, SPT-CLJ2228-5828 is a candidate for a dissociative cluster post-merger with its gas separated from dark matter, similar to the Bullet cluster (Markevitch et al. 2004). No X-ray analysis was performed until now on the system since there were no X-ray data. The cluster system was not detected either in the ROSAT All-Sky Survey (Voges et al. 1999) nor in the most recent eROSITA All-Sky Survey cluster sample (Bulbul et al. 2024). In this work, we perform the first X-ray analysis of the SPT-CLJ2228-5828 cluster system using newly obtained deep XMM-Newton data.

3. X-ray and optical data

3.1. XMM-Newton

The XMM-Newton data used in this work (Obs. ID 0884710101) were obtained during the observing cycle AO20 with total duration time of 99.1 ks (PI: H. Zohren). This was a dedicated X-ray follow-up observation to better resolve the ICM properties of the system compared to the preexisting SZ analyses.

3.1.1. Data reduction

Firstly, we reprocessed the observation data files (ODF) with the HEASoft 6.29 and XMMSAS v19.1.0 software packages. We used the emchain and epchain tasks to generate calibrated event lists. We selected only good events from the calibrated event lists, excluding bad pixels. The data were then treated for solar flare contamination. We created individual binned light curves in the 0.3−10 keV band from the entire field of view (FOV) for all three EPIC cameras (MOS1/MOS2/PN). We used 52 s and 26 s bins for MOS and PN cameras respectively. We discarded all time intervals exceeding a ±3σ Poissonian threshold of the counts bin distribution with respect to the mean. The clean exposure time of the observation is 85.8 ks for MOS1 and MOS2 and 80.1 ks for PN.

To address any residual soft proton contamination, we estimated the IN/OUT ratios (De Luca & Molendi 2004) for each camera, following Migkas et al. (2024). In a nutshell, the IN/OUT ratio refers to the surface brightness (SB) within the FOV after masking the central 10′ over the SB of the unexposed corners. Both SB values were measured in the 6 − 12 keV band for the MOS detectors and the combined (5 − 7)+(10 − 14) keV band for the PN detector. Given the decrease of the effective area beyond 10′ at > 5 keV energies, the FOV signal is not expected to have significant contributions from sky photons (e.g., cluster photons). The IN/OUT ratio was estimated using both observation data (IN/OUT_obs) and filter wheel closed (FWC) data (IN/OUT_FWC). By dividing IN/OUT_obs with IN/OUT_FWC we obtained the final IN/OUT_final ratio value. MOS and PN detectors with IN/OUT_final < 1.15 are generally considered sufficiently filtered from solar flares. For this work’s observation, all three EPIC detectors showed IN/OUT_final < 1.04, which revealed a complete lack of residual solar flare contamination.

Following the methodology of Ramos-Ceja et al. (2019) and utilizing the count-rate ratios between the 2.5 − 5 keV and 0.5 − 0.7 keV bands and their linear relation, we also ensured that our X-ray data do not suffer from Solar Wind Charge Exchange contamination (e.g., Carter & Sembay 2008). We then tested the so-called anomalous state of the MOS CCDs (Kuntz & Snowden 2008), which revealed that CCD no. 5 of MOS1 is affected by the low energy enhancement associated with the anomalous state, and thus it was discarded from the rest of the analysis. Out-of-time events were obtained and added to the instrumental background (more details below) and eventually subtracted from the data. The masking of point sources and extended sources not related to the cluster merger took place in two steps; firstly, an automatic, wavelet-based source detection algorithm is applied to the co-added MOS1/MOS2/PN (EPIC) image in the 0.5−2 keV band and masked the detected sources of no interest, as described in Pacaud et al. (2006). Then a visual inspection of the 0.5−2 keV and 2−10 keV images followed where we manually masked any suspected point sources missed by the wavelet algorithm. The masked regions were then removed from the event lists and were not considered for the spectral analysis.

3.1.2. Instrumental background

XMM-Newton data show a non-negligible level of instrumental, particle-induced background (PIB) that is important to be accounted for properly, especially when low SB systems are studied. The PIB consists of a non-vignetted continuum and fluorescence lines that originate from the interaction between high-energy particles and telescope detectors.

To estimate the PIB of the XMM-Newton data, we utilized FWC observations that are dominated by the PIB. The FWC data were rescaled to match the count-rates of the individual unexposed EPIC corners (> 925″ from the pointing’s center) of the actual observation. The energy bands used for the rescaling of the FWC data were 2.5 − 5 keV for both MOS and PN cameras and 8 − 9 keV additionally for the MOS cameras. These energy bands are free from fluorescence lines, allowing for an easier estimation of the necessary rescaling factor. For each MOS corner and PN quadrant, the PIB was taken from their closest unexposed corners. This procedure allowed us to generate PIB images that can be subtracted from the total count images, eventually resulting in PIB-subtracted, count-rate images of the source. Furthermore, PIB spectra for each EPIC detector are produced and utilized in the spectral analysis later on.

Finally, the inner unexposed PN corners show some residual contamination from FOV photons and soft protons, as demonstrated in Marelli et al. (2021). This can lead to an overestimation of the PN PIB and an underestimation of the spectroscopically measured X-ray temperature (T_X). Rossetti et al. (2024) showed that the 10 − 14 keV count-rate at > 905″ from the FOV’s center can be used to accurately rescale the FWC data and recover the unbiased PIB of PN. Additionally, Migkas et al. (2024) showed that for PN data with low IN/OUT_final < 1.1, using the 2.5 − 5 keV band at > 925″ from the FOV’s center to determine the rescaling factor of the FWC data and the PIB, also returns unbiased results for PN. Thus, the much lower IN/OUT_final values of the data used in this work guarantees a robust estimation of the PN’s PIB.

3.1.3. Imaging analysis and cluster profiles

To obtain clean count-rate images of the cluster merger we followed all the commonly used steps. Briefly, we firstly subtracted the combined EPIC PIB image from the respective combined source counts image. This resulted in ≈3900 EPIC clean source counts for each cluster in the 0.5 − 2 keV band. Then, we divided the produced image with the effective exposure map, which accounts for the absolute exposure time and the vignetting of the detectors due to the decrease of the effective area with the off-axis angle. The effective exposure map additionally accounts for the different effective areas of the PN and MOS cameras for each given energy band. This results to the final, absorbed, combined EPIC count-rate image. To acquire the respective unabsorbed image, we correct the former for the X-ray Galactic absorption using the afore-mentioned N_tot value and the energy range of interest. We check if N_tot significantly varies within the XMM-Newton FOV but we do not find any significant changes. The final, fully corrected, smoothed X-ray image of the cluster system in the 0.5 − 2 keV band is shown in Fig. 1.

Fig. 1. Fully corrected, background-subtracted XMM-Newton smoothed EPIC count-rate image of the SPT-CLJ2228-5828 cluster pre-merger in the 0.5 − 2 keV band. The R500 and R200 radii of each cluster are displayed with white dashed and cyan circles respectively. The region between the solid green circles represents the area of the sky background estimation. The emission limit of the nearby active galactic nucleus and its host galaxy is also marked with a white dashed circle.

To extract and fit the SB and gas density (which is meant to be the electron number density n_e throughout the paper if not stated otherwise) profiles, we used the pyproffit package (Eckert et al. 2020), which properly deconvolves the X-ray cluster profiles for the limited point spread function (PSF) of XMM-Newton. By integrating the SB profile over the full cluster area, one computes the X-ray flux (f_X) and luminosity of a cluster (L_X). In addition, pyproffit applies a multi-scale decomposition of the profiles to compute the 3D deprojected n_e profile. By integrating the latter profile over the full cluster volume, one obtains the total cluster gas mass M_gas.

To properly study the gas bridge between the two clusters and distinguish it from projection effects due to the superposition of the two cluster outskirts (and signal blending due to the PSF of XMM-Newton) we utilized simulated images of the system. Specifically, we firstly simulated two clusters at the positions of the real clusters by adopting the deconvolved SB profiles obtained by pyproffit in the 0.5 − 2 kev band. We simulated the cluster emission out to 1.5 R₂₀₀. Moreover, we added a point source (which represents the bright active galactic nuclei, AGN, observed between the two clusters, see Sect. 5.1) with the same flux as measured by the analysis of the real data. We also uniformly added the cosmic X-ray background (CXB) level across the full FOV. We did not add any gas bridge emission. Then, we convolved the simulated image with a 15.5″ PSF, which accounted for the half-energy width of XMM-Newton. Given that the entire cluster system is located close to the FOV center, there are no significant variations of the XMM-Newton PSF across the cluster merger. This PSF-convolved simulated image (bottom panel of Fig. 7) illustrates what one would observe with XMM-Newton purely due to the emission of the two clusters and their projection effects, assuming no gas bridge. Comparing this simulated image to the real one, we can deduce the existence of excess emission from the assumed gas bridge and infer the bridge’s X-ray properties. Finally, to check the robustness of our method, we measured the PSF-convolved SB profiles of the two clusters from the simulated image. We retrieved very similar results to the true observed profiles, which highlights the consistency of the simulated and the real image.

3.1.4. Spectral analysis

For the spectral analysis, XSPEC 12.12.0 (Arnaud 1996) was used. We consider three types of spectra: the CXB, the PIB, and the source spectra. To extract the CXB spectra, we apply two circular masks of 5′ radius centered in the X-ray peak of each cluster. This radius approximately corresponds to 2 R₂₀₀ for each cluster. We also mask the outer 12′−14′ annulus of the XMM-Newton FOV to avoid noisy photon statistics from regions with low effective area. Considering also the point source masks, we are left with ∼55% of the entire FOV to constrain the CXB, which is sufficient to provide an accurate and precise characterization of the CXB (see Fig. 1). Moreover, the PIB spectra were subtracted from the CXB and source spectra. To account for an imperfect subtraction, we added additional Gaussian lines to the model, representing Al K-α and Si-K-α fluorescence line residuals at 1.485 keV and 1.7397 keV respectively. We use the 0.5 − 7 keV energy band for the spectral analysis, avoiding the additional PN fluorescence lines in the 7.5 − 9 keV band (the high-z of the clusters, and moderate temperature as shown later, also mean there are nearly no detected source photons at ≳7 keV). The full spectral model is

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathtt{Model}=\mathtt{constant}\times (\mathtt{ape}{\mathtt{c}}_1+\mathtt{tbabs}\times (\mathtt{ape}{\mathtt{c}}_2+\\ \hfill \mathtt{powerlaw})+\mathtt{tbabs}\times \mathrm{ape}{\mathtt{c}}_3).\end{array}\end{array}$$(1)

The constant (arcmin²) term represents the scaling of the CXB components to the areas of the source regions. The components consist of the unabsorbed thermal emission from the Local Hot Bubble (apec₁ ) with a temperature free to vary within T_X ∼ 0.08 − 0.12 keV² and metallicity Z ∼ 0.9 − 1.1 Z_⊙, the absorbed Milky Way Halo emission (apec₂) with T ∼ 0.18 − 0.27 keV and Z ∼ 0.9 − 1.1 Z_⊙ (McCammon et al. 2002), and the cosmic X-ray background from the unresolved point sources (powerlaw) with a photon index of 1.46 (e.g., Luo et al. 2017). The Galactic absorption along the line of sight is modeled by a tphabs model and its total hydrogen column density parameter is set to N_H, tot = 1.89 × 10²⁰ cm⁻² according to Willingale et al. (2013). The last term accounts for the source spectra and it is an absorbed thermal emission component (tbabs×apec₃). All properties of each model component (e.g., normalizations, temperatures, and metallicities) were linked across the three EPIC detectors.

Firstly, we fit only the CXB spectra. Then, we use the best-fit values as priors for the CXB model components and we leave the full model free to vary (apec1, apec2, powerlaw and apec3) while fitting the source spectra that are composed of the source+CXB photons. Cash statistics (cstat in XSPEC) were used for constraining the best-fit parameters for the model.

In this work, we extract T_X profiles of the two clusters and the gas bridge with a typical annulus or box width of ≈30 − 32″. Owing to the finite PSF of XMM-Newton, photons originating from one region might be falsely detected in another, nearby region. This might lead to an overestimation of T_X for outer, cooler regions to which photons from hotter regions leaked. To account for this effect, we utilize crosstalk ancillary response files (ARF) between different annuli, created by the SAS task arfgen. These crosstalk ARF account for the contribution of each region to the emission of all the other regions and vice versa. All regions are then fitted simultaneously. This correction effectively removes the photon leakage between different regions due to the limited XMM-Newton PSF.

As a test, we also performed the spectral analysis for each individual cluster annulus individually without applying the crosstalk correction. The T_X changes are minimal and much smaller than the statistical uncertainties (typically ≲0.3σ shifts were found). Thus, even though the crosstalk correction improves the robustness of our analysis, it is not crucial and does not change the scientific conclusions. This was also confirmed by Bartalucci et al. (2017), who used even narrower cluster annuli than in this work.

3.2. Hubble Space Telescope observations

3.2.1. Data description and reduction

The SPT-CLJ2228-5828 system was observed by the HST as part of program GO 16488 (PI: Hannah Zohren) using ACS/WFC on March 27th and 28th, 2021, followed by observations obtained with WFC3/IR on March 31st, 2021. All observations were conducted during continuous viewing zone (CVZ) opportunities, which allow for longer total integration times per orbit. The ACS observations form a 2 × 2 mosaic, with eight dithered exposures obtained per pointing in the F606W and F814W filters each, resulting in added integration times of 3.57 ks per pointing and filter. During parts of CVZ orbits HST is observing close to Earth’s illuminated limb, leading to increased background levels due to Earthshine. We preferentially scheduled F814W observations during these orbital phases with increased background, while F606W observations are preferentially obtained during orbital phases with low background. This maximizes the depth of the F606W images, which are used for WL shape measurements (see Sect. 3.2.2). In contrast, the F814W observations are only incorporated into the photometric source selection, which has weaker depth requirements.

The WFC3/IR images were taken in the F140W filter. They consist of 12 partially overlapping exposures that cover the central regions of both cluster components, with a combined integration time of 3.79 ks.

In order to have a homogeneous depth across for the mosaic image area for the WL analysis, we do not incorporate shallower ACS F606W Snapshot imaging into our analysis, which exist for the central cluster region from program SNAP 13412 (see Schrabback et al. 2021).

The data reduction closely follows Zohren et al. (2022), which we refer the reader to for details. In summary, the reduction is largely based on the standard ACS (CALACS³) and WFC3 (calwf3⁴) calibration pipelines, which we complement with the algorithm for the correction of charge-transfer inefficiency in the ACS data described by Massey et al. (2014), as well as scripts for refined masking from Schrabback et al. (2010). The image alignment and stacking makes use of the DrizzlePac⁵, as described by Zohren et al. (2022).

3.2.2. Weak lensing measurements

Using the ACS F606W (V) images we detect objects using SExtractor (Bertin & Arnouts 1996) and measure WL galaxy shapes using the KSB+ formalism (Kaiser et al. 1995; Luppino & Kaiser 1997; Hoekstra et al. 1998) employing the implementation from Erben et al. (2001) and Schrabback et al. (2007). We model the spatially and temporally varying ACS point-spread function employing the principal component analysis model calibrated by Schrabback et al. (2010, 2018) using stellar fields. The employed KSB+ pipeline was calibrated on custom ACS-like image simulations by Hernández-Martín et al. (2020), including a correction for multiplicative shear measurement bias that depends on the galaxy signal-to-noise ratio. For the source selection we incorporate the ACS F814W (I) images, preferentially selecting background galaxies via magnitude dependent cuts in V − I color, as detailed in Schrabback et al. (2018, 2021). Removing faint sources by requiring V₆₀₆ < 26.5 AB mag and S/N_flux > 10 (define as the ratio of the FLUX_AUTO and FLUXERR_AUTO parameters from SExtractor), and applying both shape measurements cuts and the color selection, our final WL source catalog comprises 1047 galaxies within the ∼6.4′×6.4′ area covered by the ACS mosaic, corresponding to a density of 25.6 galaxies/arcmin². We do not have enough photometric bands to infer reliable photometric redshifts for individual sources in the field. Instead, we follow standard procedures to infer the source redshift distribution from well-studied external references fields, to which statistically consistent source selections are applied. In particular, our estimation of the source redshift distribution follows Schrabback et al. (2021) and is based on deep data from the CANDELS fields (Grogin et al. 2011), employing the recalibrated photometric redshift catalog from Raihan et al. (2020) that makes use of photometric measurements from Skelton et al. (2014). This yields a magnitude-dependent estimate for ⟨β⟩ and ⟨β²⟩, where the lensing efficiency β is defined through $\beta =\frac{D_{\mathrm{ls}}}{D_{\mathrm{s}}}H(z_{\mathrm{s}}-z_{\mathrm{l}})$. Here, D_s and D_ls are the angular diameter distances between the observer and the source, and the lens and the source, respectively. The Heaviside step function H(x) is unity for x > 0 and zero otherwise. Weighted according to the magnitude-dependent WL shape weight w their estimated mean values are ⟨β⟩ = 0.420 and ⟨β²⟩ = 0.205 (when assuming our reference cosmological model).

3.2.3. Mass reconstruction

The WL shear γ and the convergence κ are both second-order derivatives of the lensing potential (e.g., Bartelmann & Schneider 2001). This relationship allows for a reconstruction of the κ distribution from the shear field up to a constant offset, which is often referred to as the mass-sheet degeneracy (e.g., Kaiser & Squires 1993; Schneider & Seitz 1995). We conduct such a reconstruction based on our WL catalog, fixing the mass-sheet degeneracy by assuming a mean convergence of zero within the field and employing the algorithm from McInnes et al. (2009) and Simon et al. (2009), which incorporates a Wiener filter. In order to estimate a signal-to-noise ratio (S/N) map of the κ reconstruction, we compute the root-mean-square (rms) image of 500 noise reconstructions that are based on randomizing the ellipticity phases. The resulting S/N reconstruction is shown in Fig. 2 overlaid on the HST color image as white contours. In the image, the western component of the cluster is detected with peak significance of S/N_peak, W ≃ 4.8, whereas the eastern component is only tentatively detected with S/N_peak, E ≃ 2.3. In the reconstruction, both peaks show a tentative connection at low significance (S/N ∼ 1.5 − 2.0). We also note an extension of the western component to the northeast, coinciding with the location of several potential cluster member galaxies.

Fig. 2. Overlay of the signal-to-noise ratio map of the WL mass reconstruction (shown as white contours starting at level S/N = 1.5 in steps of ΔS/N = 0.5) and a color image of the inner cluster regions, employing the HST F140W, F814W, and F606W images for the R, G, B channels. For the non-parametric mass reconstruction shown in this image, the positions of the WL centers were left free to vary. In contrast, for the mass constraints presented in Sect. 3.2.4 we fix them to the X-ray cluster centers. The two positions, as well as the BCG positions, agree well within ≲1σ. The cyan circle indicates the SZ centre from Bleem et al. (2015), while the magenta star marks the BCG of the western cluster component. The red contours show the smoothed X-ray MOS count-rate in steps of 1.385 × 10−5 cts.

3.2.4. Weak lensing mass estimation

The system under study is composed of two interacting clusters of galaxies. To disentangle the masses of the two sub-components from the WL signal, we simultaneously fit for two spherically symmetric Navarro–Frenk–White (NFW) models. The free parameters are the two masses (M_200, west and M_200, east) enclosed within each corresponding R₂₀₀. These masses are individually tied to the NFW concentration parameter c₂₀₀ (see below) through the use of the redshift-dependent concentration–mass relation of Diemer & Kravtsov (2015), with the corrected parameters of Diemer & Joyce (2019). The centers of the two NFW models were fixed to the coordinates of the X-ray emission peaks⁶.

The NFW model parameterizes the mass density ρ(r) at physical radius r by

$$\begin{array}{c}\hfill \rho (r)=\frac{M_{\mathrm{\Delta }}}{4\pi f(c_{\mathrm{\Delta }})}\frac{1}{r{(r+\frac{r_{\mathrm{\Delta }}}{c_{\mathrm{\Delta }}})}^2},\end{array}$$(2)

where c_Δ is the concentration parameter and f(c_Δ) = log(1 + c_Δ)−c_Δ/(1 + c_Δ). Projecting the three-dimensional density model onto the plane of the sky yields the mass surface density:

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)=2{\displaystyle {\int }_0^{\infty }}\rho \left(\sqrt{R^2+{\zeta }^2}\right)\mathrm{d}\zeta ,\end{array}$$(3)

where R is a projected radius and the integration variable ζ is in the direction of the line of sight. Analytic expressions for the projected surface density Σ, the convergence κ, and the shear γ of the NFW profile have been provided by Bartelmann (1996) and Wright & Brainerd (2000). To model the shear field given two masses and their associated centers, we use the projected surface density (Eq. 3) of the NFW model and its associated prediction for the convergence and the tangential shear γ_t at the projected radius of each measured ellipticity, excluding any ellipticities with projected radii R ≤ R_min = 0.2 Mpc. Note that there are two radii associated with each ellipticity (corresponding to each X-ray or BCG center coordinate), so that effectively two circular regions are excised. We convert the tangential shears of the model into γ₁ and γ₂ by the relation

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\gamma }_1& =-{\gamma }_{\mathrm{t}}cos(2\alpha )+{\gamma }_{\times }sin(2\alpha )\hfill \\ \hfill {\gamma }_2& =-{\gamma }_{\mathrm{t}}sin(2\alpha )-{\gamma }_{\times }cos(2\alpha ),\hfill \end{array}\end{array}$$(4)

where α is the position angle of a given background galaxy, measured from the respective center, and the cross shear γ_× vanishes by construction due to the spherical symmetry of the model.

We fit the measured reduced shear estimates provided by the galaxy ellipticities applying magnitude-dependent shear weights (Schrabback et al. 2018) to estimate masses using Markov chain Monte Carlo sampling. In each step of the chain, we also scale M₂₀₀ to M₅₀₀ to allow for a subsequent comparison to the scaling relation of Bulbul et al. (2019) in Sect. 4.7.

Because reduced shear is not additive, we first computed the convergence κ_∞ and shear γ_∞ at an infinite source distance for both components of the model, added the latter, and subsequently computed the total reduced shear of the model as

$$\begin{array}{c}\hfill g=(1+(\frac{⟨{\beta }_s^2⟩}{{⟨{\beta }_s⟩}^2}-1)⟨{\beta }_s⟩{\kappa }_{\infty })\frac{⟨{\beta }_s⟩{\gamma }_{\infty }}{1-⟨{\beta }_s⟩{\kappa }_{\infty }},\end{array}$$(5)

which accounts for the β correction of Hoekstra et al. (2000). Here, β_s = β/β_∞, with β_∞ the lensing efficiency of a source at infinite redshift. Our data are the measured ellipticities ${\widehat{g}}_{k,i}$, where k ∈ [1, 2] and the second subscript indicates the ith ellipticity. The log-likelihood that we seek to maximize is given by

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{2n}{2}ln(2\pi )-{\displaystyle \sum _{i=1}^n}(2ln({\sigma }_{{\widehat{g}}_i})+\frac{1}{2}{\displaystyle \sum _{k=1}^2}{\left(\frac{{\widehat{g}}_{k,i}-g_{k,i}}{{\sigma }_{{\widehat{g}}_i}}\right)}^2),\end{array}$$(6)

where n is the number of individual ellipticities measured and the uncertainties on the reduced shear estimates ${\sigma }_{{\widehat{g}}_i}=w_i^{-0.5}$ are computed from the shear weights w_i (Schrabback et al. 2018).

To determine a 68% credibility interval of each mass component, we take the shortest interval of the respective marginalized posterior containing 68% of the samples. The best-fit is taken at the point with the highest density in the marginalized posterior. We test the robustness of the uncertainty estimates ${\sigma }_{\widehat{g}}$ of the ellipticities by subtracting the best-fit model (in terms of g₁ and g₂) from the data, randomizing the residual ellipticities by rotation, and adding the best-fit model back in. We repeat this 100 times and determine the best-fit masses in each iteration. From the 100 pairs of best-fit masses, we computed the variance for each component. We find that the results of this procedure are consistent with the mass variances of the original fit to within 10%, indicating that the ellipticity uncertainties are robust.

3.3. Gemini

3.3.1. Observations and data processing

Gemini observations were carried out under the program ID GS-2020B-Q-131 (PI: Zenteno). The goal of this program was to characterize the dynamics of SPT merging galaxy cluster candidates at z > 0.6 reported in Zenteno et al. (2020) using the redshift information of the members’ galaxies, which includes SPT-CLJ2228-5828. We selected galaxies for spectroscopy follow-up using catalogs and images⁷ from the Dark Energy Survey Year 3 (DESY3 Sevilla-Noarbe et al. 2021). Galaxies were chosen from those located in the red cluster sequence (RCS) with (r − i)_auto ≤ ±0.22 from the best fit, and located within then known R₂₀₀ radius (≤2.85′, Zenteno et al. 2020).

To have enough cluster members and to be able to detect any deviation from a Gaussian distribution in the redshift space, about ∼70 − 80 galaxy members are required to be detected for a typical SPT cluster with M_200, c ∼ 5 × 10¹⁴ M_⊙. Therefore, we have to go as deep as M^* + 1 mag. This magnitude cut-off corresponds to an apparent magnitude of i ∼ 22 AB mag at the estimated photometric z of SPT-CLJ2228-5828 (no prior spectroscopic z was available). We included 71 galaxies with i ≲ 22 AB mag in two designed masks.

Multi-object spectroscopic (MOS) observations were performed using the Gemini Multi-Object Spectrograph (hereafter GMOS, Hook et al. 2004), mounted at the Gemini South telescope in Chile, in queue mode (program ID: GS-2020B-Q-131). The masks were observed between 23 November 2020 UT and 04 December 2020 UT, during gray time, clear skies, and with a seeing that varied between 0.7″ and 1″. The spectra were acquired with the 400 line mm⁻¹ ruling density grating (R400), binning by two in both directions, 1″ slitlets, and central wavelengths of 8600, 8700 and 8800 Åoffsets are applied to avoid loss of any line that could, by chance, end in the gaps between the GMOS CCDs). Due to the faint nature of the galaxies and the central wavelength used, the contamination by the sky lines is expected to be significant. To minimize this effect and obtain reliable redshifts, masks were observed using the nod-and-shuffle technique (Glazebrook & Bland-Hawthorn 2001) in micro shuffling mode.

All spectra and corresponding calibrations were processed using the latest version of the Gemini GMOS package inside IRAF, following the standard procedures for nod-and-shuffle observations. In summary, science spectra, spectroscopic flats, and comparison lamps are bias subtracted and trimmed. Spectroscopic flats are processed by removing the calibration units’ uneven illumination and the corresponding GMOS spectral response, leaving only the pixel-to-pixel variations. Then, the two-dimensional science spectra are flat-fielded, calibrated in wavelength, sky subtracted, and extracted to a one-dimensional format using a fixed aperture of 1.5″. With the chosen slit width and central wavelengths, the final spectra have a resolution of ∼7.3 Å, and a dispersion of ∼1.34 Å pixel⁻¹.

3.3.2. Galaxy redshifts

We estimated the redshifts of the galaxies by cross-correlating the spectra with high signal-to-noise templates using the routine FXCOR available inside the IRAF RV package. For those spectra with clear emission lines, we used the routine RVIDLINE, employing a line-by-line Gaussian fit to measure the redshifts. We were able to measure the redshift for 59 galaxies out of the 71 selected for spectroscopy (∼83% of the sample) in the field of SPT-CLJ2228-5828 (5.5 × 5.5 arcmin² GMOS FOV). There are ten galaxies for which the signal-to-noise ratio was too low to estimate the redshift correctly, and two galaxies for which the spectra were observed in both masks.

4. Analysis of the merging galaxy clusters

4.1. Visual inspection and X-ray emission peak

The clean XMM-Newton count-rate image in the 0.5 − 2 keV energy band (Fig. 1) reveals two separate clusters projected close to each other, only ≈2′ apart. Both clusters seem rather symmetrical and relaxed with obvious extended emission beyond their bright centers. From the visual inspection alone, one can determine that this is not a dissociative post-merger system as initially thought from the past SZ-optical and WL analyses, but a probable pre-merger system, or two clusters at different cosmic distances projected close to each other on the plane of the sky. Our following analysis confirms that the system is a real pre-merger and gas collision has actually started.

Furthermore, there is an apparent bridge connecting the two clusters that shows excess emission. Although this could be attributed to the superposition of the projected cluster outskirts, in Sect. 5 we confirm the existence of the excess gas bridge. Northwest of the cluster system, there is a bright, slightly extended source with an apparent size of ≈1′ (limit of emission detection compared to CXB, Fig. 1). We did not find any registered optical redshift information in the literature for this object. We performed a spectral analysis to the system finding that a pow model fits the central 20″ of the northwest source better than an apec model indicating the presence of a luminous AGN in the center. The (20″ − 50″) annulus is fitted better with an apec model with $T=0.{73}_{-0.25}^{+0.32}$ keV and $z=0.{08}_{-0.05}^{+0.09}$ with an assumed metallicity of Z = 0.3 Z_⊙. These results suggest that the source is a nearby, hot elliptical galaxy with a ∼100 kpc radius projected onto the plane of the sky⁸. Therefore, we mask it for the rest of the analysis.

We consider the cluster centers to be the X-ray emission peak of each system. To determine the emission peaks, we used the clean count-rate image (produced as described in Sect. 3.1.3) in the 0.4 − 1.25 keV energy band where the cluster signal maximizes over the CXB. We first smoothed the image by a Gaussian function with a 45″ full width half maximum (FWHM) and the initial X-ray emission peaks were determined. This intentional worsening of the resolution was adopted to ensure that the X-ray emission peaks were not selected based on a poorly removed AGN or a single bright pixel of the image. After this first peak selection, we iteratively improved the Gaussian function FWHM to look for a more precise X-ray emission peak within a circle centered at the first peak with radius five times the previous FWHM. This was done until a final resolution of 7.5″ was reached. The final X-ray cluster centers in equatorial coordinates are (RA_Eastern, Dec_Eastern) = (337.239°, −58.479°) and (RA_Western, Dec_Western) = (337.171°, −58.474°) for the eastern and western cluster respectively. Their projected angular distance is 2.14′, which corresponds to 951 kpc according to their z determined in Sect. 4.2.

4.2. Redshift determination

To determine accurate redshifts for the two clusters, we implement the spectroscopic galaxy information obtained with Gemini GMOS. For each cluster, we consider only galaxies within 1′ from the cluster’s X-ray center to avoid strong contamination from the respective neighboring cluster. The spatial distributions of the galaxies and the galaxy z histograms of the two systems are displayed in Fig. 3.

Fig. 3. Spectroscopic redshift information of the two clusters. Left: Spatial distribution of the Gemini-studied galaxies with spectroscopic z. The color bar represents the z value of each galaxy. The light green and purple circles represent the 1′ area around the eastern and western clusters (X-ray peaks), respectively, and display the regions within which galaxies are assumed to belong to the respective cluster. Galaxies with bright yellow or deep blue colors do not belong to the cluster system. Middle: Spectroscopic galaxy z histogram for the eastern cluster. The black line represents the X-ray z normalized likelihood. Right: Same as in middle panel but for the western cluster.

For the eastern cluster we find a clear peak in the galaxy distribution around z ≈ 0.75 − 0.8. Five more galaxies are projected on the same line of sight, but they are not cluster members, and therefore we excluded them. There is an additional galaxy at z = 0.734, ≈3500 − 9000 km s⁻¹ away from any other remaining galaxy of the cluster. We exclude this object since it is highly unlikely that it is a part of the eastern cluster. The assigned cluster z is the weighted mean⁹ of the remaining 10 galaxy members, which is z_east = 0.771. This is consistent with the values reported on the previous studies of the entire merging system (see Sect. 2). Finally, the cluster velocity dispersion σ_vel describes the width (standard deviation) of the galaxy velocity distribution in the cluster’s rest frame and can be used as a cluster mass proxy. For the eastern cluster galaxies we measure σ_vel = 1042 ± 290 km s⁻¹.

For the western cluster, the z distribution is slightly more ambiguous and with fewer galaxy members. Firstly, we exclude the four galaxies that do not lie within 0.7 < z < 0.84 since it is obvious that do not belong to the western cluster. A peak distribution of six galaxies remains with a weighted mean of z = 0.768. The velocity dispersion of the cluster is σ_vel = 743 ± 335 km s⁻¹.

One cannot easily determine whether the difference Δz = 0.003 between the two clusters implies a different cosmic distance or if it arises because of their relative peculiar motion. Assuming that both clusters lie at the same distance (z = 0.769), their relative peculiar velocity along the line of sight is ≈520 km s⁻¹, which is a typical peculiar velocity value. Given the orientation of the double cluster system and its pre-merging state, a significant relative velocity component might also exist in the orientation of the plane, vertical to the line of sight. This is discussed further in Sect. 5.5. On the other hand, assuming no peculiar velocities, the comoving distance between the two clusters would be ≈8.5 Mpc (considering also their angular distance). This distance would imply that no physical interaction takes place between the two systems and the apparent gas bridge is the result of projection effects. However, this scenario is strongly disfavored from the analysis presented in Sect. 5.

As a cross-check, we also constrained the X-ray redshifts (z_{X − ray}) of the two clusters. For each cluster, we simultaneously fitted the X-ray spectra obtained from their < 0.5 R₅₀₀ and 0.5 − 1 R₅₀₀¹⁰ regions with different apec models but linked z_{X − ray}. We find $z_{\mathrm{X}-\mathrm{ray}}=0.{771}_{-0.012}^{+0.016}$ and $z_{\mathrm{X}-\mathrm{ray}}=0.{756}_{-0.019}^{+0.019}$ for the eastern and western clusters, respectively, fully consistent with the optical spectroscopic z. The normalized likelihoods of the z_{X − ray} values are also displayed in the middle and right panels of Fig. 3 as black lines. The weighted mean of the best-fit z when combining z_{X − ray} and z_opt differs only by < 0.02% compared to the z_opt alone.

4.3. Radius determination

Next, we needed to estimate the R₅₀₀ radius of both clusters. To do so, we implemented an iterative process of measuring X-ray cluster properties and utilizing the cluster scaling relations presented in Bulbul et al. (2019, B19 hereafter). These authors used a high-z subsample of clusters detected in the 2500 deg² SPT-SZ survey, and as such, they studied a similar cluster population with the systems that comprise SPT-CLJ2228-5828. To estimate the R₅₀₀ of the two clusters, we used the Y_X, CE − M₅₀₀ relation of B19, where Y_X, CE = M_gas × T_X, CE. Y_X, CE is the X-ray equivalent of the Compton-y parameter and a proxy of the total gas pressure, while T_X, CE is the core-excised T_X measured within (0.15 − 1) R₅₀₀. This relation exhibits the lowest M₅₀₀ scatter among the ones studied in B19 and thus offers the best precision for constraining R₅₀₀ (since $R_{500}\propto M_{500}^{1/3}$). To measure M_gas, T_X, and Y_X, CE, we consider only the half-circle region of the clusters opposite to the merging direction and the neighboring cluster¹¹. These cluster sides appear to be relatively relaxed and not significantly affected by the ongoing merging. Firstly, we assume R₅₀₀ = 1′, compute Y_X, CE, and estimate a new R₅₀₀. We repeat this iteration until the R₅₀₀ estimation converges to < 5% from the last iteration. For the eastern galaxy cluster we find R₅₀₀ = (1.52 ± 0.05)′, or (675 ± 22 kpc), and for the western cluster we find R₅₀₀ = (1.60 ± 0.05)′, or (711 ± 22 kpc). The two merging clusters have very similar radii. Their R₅₀₀ radii overlap, but do not encompass each other’s cluster center. According to Reiprich et al. (2013), R₂₀₀ ≈ 1.538 R₅₀₀, and thus R₂₀₀ = (2.34 ± 0.08)′ (1040 ± 32 kpc) and R₂₀₀ = (2.46 ± 0.08)′ (1094 ± 32 kpc) for the eastern and western clusters, respectively. For both clusters, the R₂₀₀ radii marginally encompass the center of the other cluster.

4.4. X-ray profiles and galaxy cluster properties

Knowing the z and R₅₀₀ of each cluster, we can determine a wide range of their X-ray properties and derive their radial profiles. All quantities are measured in the non-merging half-circle region of each cluster, that is, the east half of the eastern cluster and the west half of the western cluster. We derive the SB(r), T_X(r), pressure P(r), and entropy K(r) profiles. We also determine f_X, the core-included and core-excised L_X, T_X, and Y_X, M_gas, the average metal abundance Z, the X-ray concentration c_X, the characteristic pressure P₅₀₀ and entropy K₅₀₀, and the total mass M_{500, X − ray} as determined by X-rays. All the cluster properties are shown in Table 1.

4.4.1. Surface brightness profiles

The SB profiles in the 0.5 − 2 keV band¹² were extracted up to 3.5′, which corresponds to ≈(2.2 − 2.3) R₅₀₀ ≈ (1.4 − 1.5) R₂₀₀, or ≈R₁₀₀ for both clusters (Reiprich et al. 2013). The total (cluster+CXB), PSF-convolved profiles and the deconvolved cluster-only profiles are shown in Fig. 4 for both clusters. Overall, the two systems show similar SB profiles. The eastern cluster shows a slightly more peaked central emission than the western cluster, with its emission becoming weaker than the CXB at ≈0.65 R₅₀₀. Nevertheless, some excess emission is significantly detected even at > R₂₀₀. The western cluster does not show any emission beyond R₂₀₀ and it is not as centrally peaked. However, it shows a higher SB than the eastern cluster in the (0.2 − 0.6) R₅₀₀ range while its emission also becomes lower than the CXB at ≈0.65 R₅₀₀.

Fig. 4. Surface brightness profiles (PIB-subtracted) for the eastern and western clusters (left and right panels, respectively) in the 0.5 − 2 keV band. The total (cluster+CXB) PSF-convolved profiles are displayed in black, while the cluster-only PSF-deconvolved profiles are shown in the filled blue band, which represents the 68.3% error band. The R500 and R200 values are displayed with the dashed green and blue vertical lines respectively. The CXB level is displayed with the black horizontal line.

From the deconvolved SB profiles we determine the 0.5 − 2 keV core-included luminosity within R₅₀₀ to be L_CI, east = (1.11 ± 0.05)×10⁴⁴ erg s⁻¹ and L_CI, west = (1.07 ± 0.05)×10⁴⁴ erg s⁻¹ for the eastern and western cluster respectively. Both values correspond to f_X ≈ 4.5 × 10⁻¹⁴ erg s⁻¹ cm⁻². Similarly, when only the (0.15 − 1) R₅₀₀ area is considered, we measure L_CE, east = (0.67 ± 0.05)×10⁴⁴ erg s⁻¹ and L_CE, west = (0.71 ± 0.05)×10⁴⁴ erg s⁻¹ for the eastern and western cluster respectively. The two clusters have very similar X-ray luminosities.

Integrating the deconvolved SB profiles out to 0.1 R₅₀₀ we also measure the X-ray concentration $c_{\mathrm{X}}=\frac{f_{\mathrm{X}(<0.1R_{500})}}{f_{\mathrm{X}(<R_{500})}}$, which is defined as the ratio between the emission coming from the cluster’s center, 0.1 R₅₀₀, over the total emission within R₅₀₀. For the two clusters, we find c_X, east = 0.38 and c_X, west = 0.31. According to Lovisari et al. (2017), the c_X values indicate two relaxed clusters, which is surprising given their actual merging state as shown in Sect. 5. This suggests that, in cluster pairs, the opposite sides than the ones interacting probably remain unaffected, although our sample is too small to derive robust conclusions.

4.4.2. Gas density profiles

We extracted the two n_e profiles as described in Sect. 3.1.3 and fit them using a double-β model as in Shitanishi et al. (2018):

$$\begin{array}{c}\hfill n_{\mathrm{e}}(r)={(n_{0,1}^2{(1+\frac{r^2}{r_{c,1}^2})}^{-3{\beta }_1}+n_{0,2}^2{(1+\frac{r^2}{r_{c,2}^2})}^{-3{\beta }_2})}^{\frac{1}{2}},\end{array}$$(7)

where n₀, r_c, and β, are respectively the central gas density, core scale, and slope of the two profile components. Although the double-β model in Eq. (7) is a simplified version of the n_e profiles presented in Vikhlinin et al. (2006), it fits the n_e sufficiently well.

The n_e profiles of the two clusters are displayed in the top left panel of Fig. 5. Their behavior compared to each other naturally resembles the SB profiles. Both clusters show a moderate central gas density of n_0, 1 ≈ (1 − 1.25)×10⁻² cm⁻³ with the eastern cluster exhibiting a slightly denser core, while being marginally less dense at (0.2 − 0.65) R₅₀₀ than the western cluster. The central n_e values further suggest that the two clusters are most likely relaxed according to the relaxed-disturbed classification of Lovisari et al. (2017) based on several X-ray cluster properties. The scale radius of the core component is r_c, 1 ≈ 0.145 R₅₀₀ for both systems with the scale of the second component being r_c, 2 ≈ 0.58 R₅₀₀. Despite this similarity, the second components of the n_e profiles generally differ between the two clusters, with the western cluster showing a 72% higher n_0, 2 = (1.47 ± 0.16)×10⁻⁴ cm⁻³ compared to the eastern cluster, as well as a 48% higher β₂ = 0.958 ± 0.112. The onset of the dominance of the second n_e components is obvious for both clusters as the slopes change at ≈(20 − 30)″. This further indicates the presence of relatively relaxed cluster cores. Furthermore, the density drops to n_e ≲ 2 × 10⁻⁴ cm⁻³ at > R₅₀₀ and n_e ≲ 10⁻⁴ cm⁻³ at the very outskirts (> R₂₀₀), with both clusters showing typical gas distributions. Integrating the n_e profiles, we measure the gas masses within R₅₀₀ to be M_gas, east = (1.76 ± 0.13)×10¹³ M_⊙ and M_gas, west = (2.06 ± 0.15)×10¹³ M_⊙ for the eastern and western cluster respectively. This further confirms the similarity of the two clusters and their average mass.

Fig. 5. X-ray radial profiles for the eastern and western clusters (green and blue colors, respectively). The R500 values are displayed with the dashed vertical lines. Upper left: Gas (electron number) density ne profiles with their best-fit models. The R200 values are also displayed in this panel. Upper right: Temperature profiles (projected, 2D) with their best-fit deprojected, 3D models. Due to the limited data, no model uncertainties could be derived. Bottom left: Inferred pressure profile best-fit models computed as P = ne × TX, 3D. The 1σ error band is underestimated since it accounts only for the uncertainty of the ne profile assuming no uncertainty of the TX, 3D profile. Bottom right: Inferred entropy profile best-fit models computed as $K=T_{\mathrm{X},3\mathrm{D}}\times n_{\mathrm{e}}^{-\frac{2}{3}}$. The 1σ error band is underestimated since it accounts only for the uncertainty of the ne profile assuming no uncertainty of the TX, 3D profile.

4.4.3. Temperature profiles and metal abundance

Firstly, we determine the overall T_X values for both clusters. For the eastern cluster, we find $T_{\mathrm{X},\mathrm{CI}}=3.{64}_{-0.35}^{+0.43}$ keV and $T_{\mathrm{X},\mathrm{CE}}=3.{22}_{-0.50}^{+0.59}$ keV for the core-included and core-excised regions respectively. Similarly, for the western cluster we find $T_{\mathrm{X},\mathrm{CI}}=3.{84}_{-0.35}^{+0.36}$ keV and $T_{\mathrm{X},\mathrm{CE}}=3.{41}_{-0.45}^{+0.54}$ keV. The two clusters show very similar T_X, with the western cluster being slightly hotter, which is consistent with the marginally higher values of L_CE and M_gas it also shows compared to the eastern cluster.

The universal metal abundance Z could be only loosely constrained with $Z_{\mathrm{east}}=0.{84}_{-0.37}^{+0.42}{Z}_{\odot }$ and $Z_{\mathrm{west}}=0.{76}_{-0.31}^{+0.36}{Z}_{\odot }$. Although the best-fit values are relatively high, the large uncertainties prevent us from drawing firm conclusions. For all following T_X constraints, we let Z free to vary within the 1σ limits of the < R₅₀₀ constraints.

Due to the small apparent size of the two clusters, extracting and fitting T_X profiles is challenging. To have a sufficient number of counts per bin, avoid strong photon mixing between bins due to the limited PSF, and have enough bins to constrain the T_X profile, we define three bins for each cluster: 0 − 0.3 R₅₀₀, 0.3 − 0.6 R₅₀₀, and 0.6 − 1 R₅₀₀. For each bin, we measure the 2D projected T_X. Once again, the two cluster systems show very similar behavior. A gradual T_X decrease is observed with the central, middle, and outer bins showing ≈3.8 − 3.9 keV, ≈3 − 3.3 keV and ≈2 − 2.1 keV respectively, for both clusters.

To determine the 3D deprojected T_X, 3D(r) profiles, we project the profiles onto the line of sight and fit the measured T_X. To do so, we adopt the T_X, 3D(r) profile described in Vikhlinin et al. (2006) and the projection weighting scheme introduced by Vikhlinin (2006). The T_X, 3D(r) parameters are then determined by considering the resulted projected T_X, 2D(r) profile that describes the data best. This methodology has been implemented by several past studies (see, e.g., Bartalucci et al. 2017; Schellenberger & Reiprich 2017; Lovisari et al. 2020, for more details). The T_X, 3D(r) profile is given by

$$\begin{array}{c}\hfill T_{\mathrm{X},3\mathrm{D}}(r)=T_{\mathrm{X},\mathrm{CE}}^{\ast }{T}_{\max }\frac{{(r/r_{\mathrm{cool}})}^{a_{\mathrm{cool}}}+T_0/T_{\max }}{{(r/r_{\mathrm{cool}})}^{a_{\mathrm{cool}}}+1}\frac{{(r/r_t)}^{-a}}{{(1+{(r/r_t)}^b)}^{c/b}},\end{array}$$(8)

where T_X, CE^* is the T_X measured within (0.15 − 0.75) R₅₀₀¹³. Moreover, T₀, r_cool, and a_cool drive the behavior of the T_X, 3D(r) profile close to the core (≲[0.05 − 0.15] R₅₀₀), while T_max, r_t, b, and c describe the profile further out until the outskirts. Since we only have three measured T_X bins, this eight-parameter model will clearly overfit the data. To this end, we fix (r_cool, a, a_cool, r_t, b) = (0.039 R₅₀₀, 0.036, 0.898, 1.095 R₅₀₀, 2.36) to the best-fit values obtained by Chen et al. (2024) for the EXCPReS (Evolution of X-ray galaxy Cluster Properties in a Representative Sample) catalog, which contains 0.4 < z < 0.6 clusters with similar properties as the ones we focus in this work. Also, Chen et al. (2024) constrained T₀ = 0.509 ± 0.452 but we fix it to T₀ = 0.9, which is more suitable for our two clusters¹⁴. Thus, we only fit the T_max and c parameters, which have the most impact on our T_X, 3D profiles. The projected best-fit T_X profiles, alongside the T_X measurements, are displayed in the top-right panel of Fig. 5. The statistical uncertainties of the T_X, 3D(r) profiles cannot be robustly constrained due to the limited data and the wide range of model parameters, many of which we kept fixed, which underestimates the uncertainties. Therefore, we only focus on the best-fit T_X profile functions. The two clusters return very similar results near their center (< 0.3 R₅₀₀) and toward the outskirts (> 0.65 R₅₀₀). In the in-between cluster regions the western cluster seems to be ≈11% hotter, although this is well within the statistical uncertainties. Even if very limited information on the T_X-profiles is available, the generally used profile in Eq. (8) gives a very good description of the data even with most parameters fixed to literature values. As such, to derive the pressure and entropy profiles in the next sections, we use the best-fit functions of the T_X profiles, assuming no uncertainties.

4.4.4. Pressure profiles, X-ray Compton equivalent Y_X, and total mass M_500, YX

As explained in Sect. 4.3, the X-ray equivalent for the Compton parameter is defined as Y_X = M_gas × T_X and is a proxy of the total gas pressure of the cluster. The core-included and core-excluded Y_X values are computed by the use of T_X, CI and T_X, CE respectively. For the eastern cluster, we find $Y_{\mathrm{X},\mathrm{CI}}=6.{40}_{-0.80}^{+0.92}\times {10}^{13}{M}_{\odot }$ keV and $Y_{\mathrm{X},\mathrm{CE}}=5.{67}_{-1.01}^{+1.12}\times {10}^{13}{M}_{\odot }$ keV while for the western cluster we find $Y_{\mathrm{X},\mathrm{CI}}=7.{91}_{-0.93}^{+0.96}\times {10}^{13}{M}_{\odot }$ keV and $Y_{\mathrm{X},\mathrm{CE}}=7.{01}_{-1.07}^{+1.23}\times {10}^{13}{M}_{\odot }$ keV. The total mass within R₅₀₀ as inferred from the Y_X, CE − M₅₀₀ relation of B19 (after correcting for the different cosmology) is M_500, YX = (2.10 ± 0.20)×10¹⁴ M_⊙ and M_500, YX = (2.41 ± 0.22)×10¹⁴ M_⊙ for the eastern and western clusters respectively¹⁵.

The cluster gas pressure is the least affected thermodynamic property by the thermodynamic evolution of clusters. The gas pressure profile, assuming an ideal and monoatomic gas, is given by P(r) = n_e(r)T_X, 3D(r). Since we have derived both the n_e(r) and T_X, 3D(r) profiles, we take their product to be the P(r) profiles. The latter are shown in the bottom-left panel of Fig. 5. Due to the limited T_X information, we cannot infer the P(r) profile much further than R₅₀₀. Nevertheless, we see that both P(r) profiles show a typical behavior for clusters of such T_X and z (e.g., Ghirardini et al. 2017; He et al. 2021). According to the results of Arnaud et al. (2010), the P(r) profiles indicate averagely relaxed clusters without strong cool core presence (but possible the presence of weak cool cores). Moreover, under the assumption of self-similarity, the characteristic pressure P₅₀₀ of a cluster is given by Voit (2005), Nagai et al. (2007)

$$\begin{array}{c}\hfill P_{500}=1.89\times {10}^{-3}E{(z)}^{\frac{8}{3}}\times {\left(\frac{M_{500}}{3\times {10}^{14}{M}_{\odot }}\right)}^{\frac{2}{3}}{h}^{8/3}{\mathrm{keV}/\mathrm{cm}}^3,\end{array}$$(9)

where E(z) = $\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$ and h = H₀/(70 km s⁻¹ Mpc⁻¹). For the eastern and western cluster we find P₅₀₀ = (4.92 ± 0.33)×10⁻³ keV cm⁻³ and P₅₀₀ = (5.42 ± 0.32)×10⁻³ keV cm⁻³ respectively, which approximately corresponds to the P(r) value at r ≈ 0.3 R₅₀₀, consistent with Ghirardini et al. (2017) within the scatter of clusters at similar z. Finally, we do not have sufficient data to probe the previously detected flattening of the profile at < 0.1 R₅₀₀ (e.g., McDonald et al. 2014).

4.4.5. Entropy profiles

The entropy K of the ICM is an essential tool for studying the thermodynamic state and history of the cluster gas. It always increases when gas heating occurs (e.g., Voit 2005) and can help us understand the relaxation state of the cluster core and outskirts, as well as their gas clumpiness (e.g., Bulbul et al. 2016). The entropy profile is given by K(r) = T_X, 3D(r)/n_e(r)^2/3, and thus, it can be directly inferred from the known T_X, 3D(r) and n_e(r) profiles. The K(r) profiles are shown in the bottom right panel of Fig. 5. As discussed before, they are only loosely constrained due to the limited T_X information. Nevertheless, one sees that K(r) decrease toward the cluster center at a level expected for relatively average clusters in terms of relaxation state (K ≈ 100 − 130 keV cm² at the core, e.g., Cavagnolo et al. 2009; Ghirardini et al. 2017). Toward the outskirts, one sees that the slope of the K(r) profiles flatten, which indicates that the cluster outskirts are not yet in dynamical equilibrium (although no robust conclusions can be drawn due to the large uncertainties). Under the self-similar model assumption, the characteristic entropy K₅₀₀ of the clusters is Voit et al. (2005)

$$\begin{array}{c}\hfill K_{500}=761\times E{(z)}^{-\frac{2}{3}}\times {\left(\frac{M_{500}}{3\times {10}^{14}{M}_{\odot }}\right)}^{\frac{2}{3}}{h}^{\frac{2}{3}}{\mathrm{keV}\mathrm{cm}}^2.\end{array}$$(10)

For the eastern and western cluster, we find K₅₀₀ = 451 ± 31 keV cm⁻³ and K₅₀₀ = 498 ± 30 keV cm⁻³, respectively, which approximately corresponds to the K(r) value at r ≈ (0.75 − 0.8) R₅₀₀.

4.5. Weak lensing masses

The two-component NFW fit to the WL data resulted in the following masses for the cluster merging system: $M_{200}=3.{26}_{-1.24}^{+1.25}\times {10}^{14}{M}_{\odot }$ for the western cluster, $M_{200}=0.{56}_{-0.49}^{+0.80}\times {10}^{14}{M}_{\odot }$ for the eastern cluster, and $M_{\mathrm{tot}}=3.{82}_{-0.97}^{+1.90}\times {10}^{14}{M}_{\odot }$ for the total mass of the system. The total mass of the system was determined by adding the M₂₀₀ of the cluster components for each step in the MCMC. We also obtain $M_{500}=2.{21}_{-0.85}^{+0.86}\times {10}^{14}{M}_{\odot }$ for the western cluster and $M_{500}=0.{40}_{-0.34}^{+0.53}\times {10}^{14}{M}_{\odot }$ for the eastern cluster. Here, the quoted WL mass uncertainties include shape noise only (see Sect. 4.6 regarding the impact of additional noise components).

Interestingly, the eastern cluster is marginally detected by the WL data, only at a 1.14σ level. Moreover, the western cluster appears to be approximately six times more massive compared to the eastern component (although only at 1.8σ), in contrast with the similarity the two clusters exhibit in X-rays. Finally, the total mass of the system (sum of the two M₂₀₀ values) is consistent with the past SZ constraints on the total mass as presented in Sect. 2.

4.6. Comparison of X-ray and weak lensing cluster masses

The X-ray masses were estimated from the Y_X, CE − M₅₀₀ relation of B19. The M₅₀₀ values in that work were inferred from the unbiased SZ S/N ratio ζ and the ζ − M₅₀₀ scaling relation of de Haan et al. (2016). The latter calibrated the ζ − M₅₀₀ relation using masses obtained by Y_X values measured by Chandra and the Y_X − M₅₀₀ from Vikhlinin et al. (2009), where M₅₀₀ represented the hydrostatic masses. The normalization of the Y_X − M₅₀₀ relation was adjusted in de Haan et al. (2016, see their Sect. 4.2 for more details) so the resulting M₅₀₀ would be consistent with more recent WL mass estimates from Hoekstra et al. (2015). This way they accounted for any systematic offset between the WL and the original hydrostatic M₅₀₀. Given this sequence of scaling relation calibrations and assuming no significant biases were propagated through these relations, we expect that the derived M_500, YX values in this work are consistent with the WL mass estimates.

To compare the WL and X-ray cluster masses, we assume that R₅₀₀ = 1.538 R₂₀₀ (Reiprich et al. 2013), which results in M_200, YX = (3.06 ± 0.29)×10¹⁴ M_⊙ and M_200, YX = (3.51 ± 0.32)×10¹⁴ M_⊙ for the eastern and western clusters respectively. The western cluster mass estimates agree well between WL and X-ray within < 0.2σ, with $M_{200,\mathrm{WL}}=0.{93}_{-0.30}^{+0.50}{M}_{200,Y_{\mathrm{X}}}$. On the other hand, the eastern cluster is marginally detected in WL and, as a result, the two mass estimates disagree at a 2.9σ level, with $M_{200,\mathrm{WL}}=0.{18}_{-0.16}^{+0.28}{M}_{200,Y_{\mathrm{X}}}$. Even though its statistical uncertainties are large, the M_200, WL value of the eastern cluster is surprisingly low. Due to the similarity of the two clusters in their X-ray properties, if there was a true lack of WL signal from the eastern cluster it would be striking and highly interesting for a follow-up study. However, we stress that the (already large) quoted WL mass uncertainty includes shape noise only. Additional statistical uncertainties are caused by large-scale structure projections and line-of-sight variations in the source redshift distribution. For individual clusters at the mass and redshift range of the target and analyses of similar setup these jointly amount to an additional statistical uncertainty of ΔM_200c, WL ≃ 0.8 × 10¹⁴ M_⊙ (e.g. Schrabback et al. 2018), which would need to be added in quadrature to the comparable shape noise uncertainties. In the special case of two clusters located in close proximity on the sky, these noise contribution would however be correlated. This is expected to reduce the relative impact of these uncertainties on the mass ratio of the components, but makes it difficult to quantify the impact precisely. Nevertheless, it is clear that the true significance of the discrepancy between the X-ray and WL mass estimates of the western component is lower than what is estimated above, making it plausible that this could simply be the result of a statistical fluke.

4.7. Consistency with known scaling relations of single clusters

Merging clusters are often excluded from samples constructed to constrain average scaling relations (e.g., Migkas et al. 2020) since the merging process is expected to alter the scaling relation behavior of such clusters. X-ray and SZ scaling relations are mostly driven by the gravitational effects on the ICM. At the same time, the ICM behavior in merging clusters is strongly affected by formed shock fronts and general gas disruption. However, limited effort has been shown to evaluate if the outer, undisturbed side of pre-merger systems follow an average scaling relation behavior. To assess this using the two clusters studied in this work, we estimate the consistency of their X-ray and WL properties as measured from their undisturbed sides with the scaling relations presented by B19. In that work the authors used only isolated, non-merging clusters to constrain scaling relations between X-ray observables and M₅₀₀. Using M_500, WL for these comparisons, we find that the western cluster lies within < 1σ for all seven scaling relations studied in B19. For the eastern cluster, such comparisons are not feasible due to the marginal detection of the cluster in the WL data and its likely underestimated M_500, WL value.

To compare the consistency between the clusters’ X-ray properties and the B19 scaling relations, we perform the following. Based on the B19 scaling relations, we infer M₅₀₀ with all available X-ray properties (M_gas, core-included and core-excluded L_X, T_X, and Y_X) and compare the results. For example, if M₅₀₀ from Y_X and L_X agree within the statistical scatter, then the studied cluster is considered consistent with the B19 scaling relations. That is, it would lie close to the hypothetical Y_X − L_X scaling relation based on the B19’s measurements. For the western cluster, we find that M₅₀₀ estimates from all X-ray properties agree within 1.2σ (given the M₅₀₀ scatter of each relation). The most discrepant quantities are T_X, CE and M_gas, which predict a M₅₀₀ difference of ≈20 ± 17%. For the eastern cluster, we find similar results with a general M₅₀₀ agreement within 1.1σ and the most discrepant quantities L_X, CI and M_gas predicting different M₅₀₀ by ≈30 ± 27%.

Clearly, we do not have a sufficient number of clusters to draw robust conclusions while the M₅₀₀ scatter of the B19 scaling relations is rather large, covering most of the inferred M₅₀₀ differences. Nonetheless, the results of both clusters favor the case that the undisrupted side of clusters at the early merging stages can be used to measure X-ray properties that would correspond to isolated, undisturbed clusters. This potentially enables the use of such merging clusters as single, isolated clusters, as long as their merging side is ignored.

4.8. Consistency with the SPT-SZ observations and cosmological implications

As discussed in Sect. 2, the past SZ studies of the merger by Bleem et al. (2015), Bocquet et al. (2019), Hilton et al. (2021) detected a single cluster with the gas thought to be between the two X-ray peaks. Therefore, we wish to determine if the SZ image is consistent with the X-ray data revealing two merging clusters, or if there is evidence for unknown systematic biases between the SZ and X-ray datasets. To do so, we use the known X-ray properties of the clusters to construct a simulated y-image with the angular resolution of SPT-SZ at 150 GHz and compare it with the real image. Specifically, we utilize the known deprojected n_e(r) and T_X, 3D(r) to create a 3D y_SZ model of the system, where y_SZ ∝ ∫_losn_e(r)T_X, 3D(r). We consider the best-fit n_e(r) and T_X, 3D(r) profiles of the two clusters and the respective values for the gas bridge, to fully account for any possible y_SZ signal. We then integrate along the line of sight to obtain a simulated y_SZ image with perfect resolution. Finally, we convolve the image with a Gaussian kernel of FWHM = 1.1′ to reproduce the expected SPT-SZ image at 150 GHz based on the X-ray data¹⁶. The resulted simulated SPT-SZ image and the real, observed image of the system by SPTpol are displayed in Fig. 6. The two images are highly consistent. This indicates that the sole reason for the false inference of the gas position by the past studies is the large angular resolution of SZ surveys.

Fig. 6. Image from SPTpol of the cluster merger (top) and the expected SPT-SZ Compton-y image based on the derived X-ray ne and TX profiles and the SPT-SZ angular resolution (bottom). The two images are consistent. That is, the SZ data are consistent with the X-ray analysis of this work.

The SZ cluster counts are extensively used for obtaining cosmological constraints (e.g., Planck Collaboration XXIV 2016; de Haan et al. 2016; Bocquet et al. 2019, 2024). As shown in this work, SZ data might confuse individual merging clusters for single objects. This might result in underestimating the total mass of such systems. As an example, for SPT-CLJ2228-5828, the mass estimates from both the SPT-SZ (Bleem et al. 2015; Bocquet et al. 2019) and ACT (Hilton et al. 2021) surveys underestimate the total mass of the cluster merger. Specifically, SPT-SZ returns ∼50% and ∼30% lower total mass than the X-ray and WL estimates respectively, while ACT returns ∼24% and ∼11% lower masses respectively. When cosmological forward modeling is utilized, one obtains an observed cluster sample from a simulated mass-selected cluster sample folded with the selection function. The latter should account for beam size effects that might lead to mislabeling merging clusters, such as SPT-CLJ2228-5828. If cosmological simulations accurately recover the fraction of such neighboring cluster systems compared to real data and the selection function efficiently describes what one would observe in the case of such systems, no cosmological bias is expected. This needs to be thoroughly checked when the beam size of the instrument used does not allow for the disentangling of such structures. On the other hand, if the above conditions are not met, then some propagated bias to the constrained cosmological parameters might be expected. However, such systems are probably rare since when X-ray followup studies are performed (e.g., B19 Lovisari et al. 2020), they do not often reveal such missed mergers. Hence, it is not expected that clusters closely projected on the plane of the sky will significantly alter the cosmological constraints of SZ cluster studies.

5. Analysis of the gas bridge between the clusters

From the X-ray image of the system (Fig. 7), the two clusters seem to interact. However, one cannot immediately distinguish if the image is solely the result of the projection of two clusters at different distances and no interaction takes place. In this case, the apparent emission between the clusters would be caused by the superposition of their outskirt emission. To determine if there is an actual gas bridge one needs to compare the observed emission with the expected emission if all the X-ray signal was originating from the summed cluster outskirt emission. In the following sections we extract the SB, T_X, P, and K profiles of the region to understand the nature of the existing gas there.

Fig. 7. Zoomed-in real (top) and simulated (bottom) image of the system in the 0.5 − 2 keV band. The simulated image assumes only the emission of the clusters and no gas bridge. The displayed boxes represent the regions in which the SB and ne in Figs. 8 and 9 were measured.

5.1. Surface brightness profile of the bridge

Firstly, we extract the SB profile along the axis of the apparent bridge using 19 box regions, as displayed in the top panel of Fig. 7. Each box has a 90″ height, 15″ width, and an 8° tilt. There is an AGN located slightly south from the bridge’s center at (RA, Dec) = (337.200° , − 58.480°) that we mask with a 15″ radius circle (no z information is available for this object and one cannot tell if it is part of the merger). Additionally, we create a simulated image of the merging system inputting the derived SB cluster profiles (Sect. 4) and the CXB level. We assume no excess emission from the bridge on top of the cluster emission. The simulated image was convolved with the XMM-Newton EPIC PSF to achieve the same angular resolution as in the real data. We also accurately simulated the AGN emission and masked it as in the real image. The simulated image is displayed in the right panel of Fig. 7. Using the simulated image, we extract the SB for the boxes between the two clusters. These values correspond to what one would measure purely due to the sum of the clusters’ emission assuming no extra gas.

In Fig. 8, we compare the simulated SB values with the true SB across the axis of the bridge. This comparison reveals that there is significant excess emission from the center of the bridge compared to the expected SB from the cluster superposition only. Specifically, there is an SB excess of (287 ± 49)% in the central box of the bridge at a 5.8σ level. Moreover, there is a ≈(57 − 65)% SB excess in the two neighboring boxes, 15″ (111 kpc) away from the bridge’s apparent center. As we move further away from the center of the bridge and toward the cluster cores, there is ≈40% excess emission toward the eastern cluster 30″ (222 kpc) away from the center of the bridge while there is no apparent excess toward the western cluster at the same distance from the bridge’s center. These results strongly demonstrate the existence of a gas bridge between the two clusters with significant X-ray emission.

Fig. 8. Surface brightness profile (CXB-subtracted) across the merging axis for the boxes shown in Fig. 7, measured in the 0.5 − 2 keV band. The black points represent the measured SB values. The gray ones represent the expected SB values from cluster projection effects and assuming no excess gas bridge. The position of the eastern and western X-ray peaks are displayed by the dashed green and blue vertical lines respectively. The horizontal black dashed line shows the CXB level. The SB ratio between the measured values and the expected values if there was no gas bridge is shown in the bottom panel.

The formation of this bridge could follow two scenarios. Firstly, a gas filament could originally connect the two clusters, and as they approached each other due to their mutual gravitational attraction, some of the filamentary gas compressed to form the highly emitting gas bridge. If the compression is substantial, weak axial shock fronts could also arise, heating the filamentary gas, similar to Sarkar et al. (2022). This scenario would not require the true R₅₀₀ of the clusters to overlap. The second case is that the material of the two clusters has started to collide due to the cluster merging, possibly forming gas clumps and stronger shocks and heating the gas. Such shocks result in significantly higher emission compared to the cluster outskirts. In that case, the cluster material could additionally mix with the gas from the presumed filament that originally connected the two systems. The two scenarios are not easy to distinguish, given the observational limitations. Thus, we study additional X-ray properties of the bridge to further understand the exact physical process that formed it.

5.2. Temperature profile of the bridge

To extract the T_X profile across the bridge axis, we need more photons per box than for the SB profiles. Hence, we consider boxes with 31.5″ width and 90″ height, which results in three equally sized, independent T_X bins between the two cluster centers. We cannot accurately predict the expected T_X if there was no gas bridge; therefore any such comparison is approximate. The T_X profile is displayed in the top left panel of Fig. 9. Across the gas bridge, T_X shows a significant increase instead of gradually decreasing as one would expect if the X-ray signal was coming only from the cluster outskirts. Specifically, at the bridge’s center, we measure $T_{\mathrm{X}}=5.{53}_{-1.09}^{+1.47}$ keV instead of the expected T_X ≈ 1.5 − 2 keV. The metallicity cannot be constrained and results in $Z=0.{29}_{-0.29}^{+0.94}{Z}_{\odot }$. Therefore, for the bridge box bins we fix it to the typical value of Z = 0.20 Z_⊙, as derived by Lovisari & Reiprich (2019) for cluster outskirts¹⁷.

Fig. 9. X-ray profiles across the cluster merging axis for the temperature (top left), gas (electron number) density (top right), pressure (bottom left), and entropy (bottom right). The TX, P, and K profiles are measured for larger bins than the ones in Fig. 7, with 31.5″ width as explained in Sect. 5.2. The gray points in the ne profile describe the expected ne value due to projection effects if no gas bridge was present. All X-ray properties are elevated compared to the case where no gas bridge or shock front exists between the two clusters.

The position of the T_X jump coincides with the SB jump in Fig. 8 and it demonstrates that the gas there has been heated by a physical mechanism. The most likely scenario is that the hot ICM of the two clusters and the gas bridge has started colliding and compressing, giving rise to a possible shock front, ≈475 kpc away from each cluster’s core (assuming the two clusters lie at the exact same cosmic distance). Moreover, this T_X jump strongly disfavors the first scenario discussed in Sect. 5.1 where the gas bridge was purely coming from a pre-existing filament and the true radii of the clusters did not interact. Such a high T_X is not expected without the interaction of the already-hot ICM of the two systems. Thus, one concludes that the true distance of the clusters should be smaller than the sum of their R₅₀₀. Considering their angular distance, this implies that their radial distance should be ≲1 Mpc. As a result, one can safely assume that the clusters lie at the same cosmic distance and their Δz = 0.003 can be attributed to a relative 520 km s⁻¹ peculiar velocity. Finally, the bin between the bridge center and the eastern cluster also shows an elevated T_X compared to the symmetrical bin opposite of the eastern cluster’s center. This suggests that any gas interaction and shock front has likely started to extend in that region.

5.3. Gas density profile of the bridge

The n_e values for the boxes can be derived from the spectral fits for each region, namely the normalization value of the apec model that describes the thermal emission of the clusters or the bridge. However, this does not allow us to disentangle the n_e coming from the cluster outskirts and the n_e coming from the gas bridge overdensity not associated with the clusters. Consequently, we use a novel approach to estimate the extra n_e originating from the bridge only for each box. Our approach is based on the already extracted n_e profiles from Sect. 4.4.2. For the following procedure, we assume both clusters to be at the same distance (z = 0.77) as discussed in the previous section.

We create a 3D model of each cluster using the known cluster n_e profiles, assuming no gas bridge or shock front. We then find the corresponding cluster volume covered by each projected box and measure the average n_e in the box based on the 3D n_e model. We integrate across the length of each box, which we consider to be from the foreground to the background (compared to each cluster’s center) R₂₀₀ limits of each cluster. The chosen limits directly and strongly affect the estimated average n_e. The latter drops significantly toward the outskirts and, the narrower the box length selected (i.e., integration limits), the higher the final average n_e. Hence, it is important to stress that the box n_e values are partially arbitrary and they are hard to physically interpret, as they cover multiple cluster layers¹⁸. Therefore, they should not be directly compared with the n_e profiles displayed in Fig. 5 that correspond to annuli at a fixed distance. The box n_e would naturally be lower. However, they allow us to compare the expected projected n_e in the bridge region with the observed n_e and determine the gas overdensity and length of the bridge. A similar approach for studying similar pre-merger systems was applied by a few past studies (e.g., Omiya et al. 2023).

Due to the coverage from each projected box of several cluster layers and the gas bridge, it is challenging to directly estimate the observed n_e from the SB values. Thus, we follow an empirical approach. We use the SB and n_e values of the 12 independent boxes covering the areas from the cluster X-ray peaks toward their outskirts, opposite to the merging direction. We fit the relation SB_i/SB_j = A(n_e, i/n_e, j)^B where the indexes i and j indicate different box regions. To fit the relation we use all the available box permutations for which SB_i/SB_j < 5, comparable to the $\frac{\mathrm{observed}\mathrm{bridge}\mathrm{SB}}{\mathrm{predicted}\mathrm{bridge}\mathrm{SB}}=3.87$ ratio from Sect. 5.1 and avoiding more extreme SB differences. From the obtained best-fit relation, we interpolate the true n_e of the gas bridge boxes by considering the $\frac{\mathrm{observed}\mathrm{bridge}\mathrm{SB}}{\mathrm{predicted}\mathrm{bridge}\mathrm{SB}}$ and the expected n_e per box as determined earlier. We compute the statistical uncertainty of the true n_e estimates through error propagation, considering also the (low) scatter of the fitted SB–n_e relation. To account for any neglected systematics in our methodology we add in quadrature an additional 5% uncertainty in all interpolated n_e values (which is a rather conservative approach). The comparison between the observed and the expected n_e per box is shown in the top right panel of Fig. 9.

Similarly to the SB comparison, there is a gas overdensity compared to the expected n_e purely from the cluster outskirts. Specifically, the central bridge region shows a (64 ± 22)% n_e excess at a 2.9σ level¹⁹. Moreover, there is a gas density jump at the bridge’s center compared to the nearby regions, further supporting the existence of a shock front. The density of the extra gas at the center of the bridge, on top of the cluster gas, seems to be n_e = (3.93 ± 1.36)×10⁻⁴ cm⁻³. This corresponds to a baryon overdensity of δ_b ≈ 320, which hints at a denser gas distribution than typical bridges between clusters, although this could be an overestimate due to the suspected shock front. The total gas density (bridge+cluster gas) is n_e = (1.01 ± 0.15)×10⁻³ cm⁻³, which corresponds to δ_b ≈ 825. Except for the central part of the bridge, the neighboring box regions also show gas overdensities. The extra n_e at ≈111 kpc (15″) from the center of the bridge is n_e ≈ 1.06 − 1.22 × 10⁻⁴ cm⁻³, corresponding to δ_b ≈ 87 − 98. Furthermore, at 222 kpc from the bridge’s center and toward the eastern cluster, the extra n_e ≈ 9.4 × 10⁻⁵ cm⁻³ and δ_b ≈ 74. Thus, we conclude that the length of the bridge formed by material that does not come from the cluster outskirts is ≈333 kpc (45″). Finally, from the estimated n_e we can compute the M_gas of the gas bridge material only. We assume a cylinder shape for the bridge with 1040 kpc radius (equal to the eastern cluster’s R₂₀₀) and 333 kpc length (covering the four boxes for which we see a n_e overdensity). We find M_gas = (1.37 ± 0.34)×10¹² M_⊙. Finally, given the measured T_X, we can also estimate Y_X = (7.57 ± 2.3)×10¹² M_⊙ × keV.

As a cross-check, we additionally derive the total n_e (clusters+bridge) using the spectral results and the apec normalization from the central T_X bin of the bridge. We assume the same cylinder shape for the bridge as before. For converting the apec normalization to n_e we follow (Veronica et al. 2022, see their Eq. 18 for details). For the 333 kpc of the bridge we obtain a total n_e = (9.94 ± 1.48)×10⁻⁴ cm⁻³, which is fully consistent with the result from our empirical approach above and strongly highlights the robustness of our employed methodology used to derive the extra n_e coming purely from the gas bridge.

5.4. Pressure and entropy profiles of the bridge

The P(r) and K(r) profiles of the bridge are displayed in the lower panels of Fig. 9. As they are both a combination of the T_X and n_e profiles, they are rather noisy. Since we do not have a predicted T_X for the case of no gas bridge, there is also no such prediction for P and K. Nevertheless, one sees that P shows elevated values across the entire bridge and an increase at the central bridge position compared to what one would expect based on the anti-symmetric P values with respect to the cluster centers. Characteristically, 60″ away from the cluster centers and opposite to the bridge, we measure P ∼ (6 − 7)×10⁻⁴ keV cm⁻³. Assuming the same T at the same distance toward the bridge, and given that P = n_e × T_X, one would expect to measure P ∼ (1 − 1.5)×10⁻³ keV cm⁻³ at the middle of the bridge. However, we measure $P=5.{12}_{-1.12}^{+1.44}\times {10}^{-3}$ keV cm⁻³, which constitutes a substantial P increase at the position of the suspected shock front. Similarly, entropy increases at the bridge’s center with $K={583}_{-152}^{+169}$ keV cm², showing larger values than the cluster outskirts opposite to the merging direction as presented in Fig. 5, although the difference is not statistically significant due to the large uncertainties. Even if not much information can be robustly drawn from the P and K profiles due to the complex geometry of the boxed covered area, the results further support the existence of a weak-to-moderate shock front.

5.5. The shock front and the cluster collision geometry

The detection of shock fronts in cluster (pre-)mergers is relatively rare due to the many requirements that need to be met for such a shock to be formed and observed. From Fig. 9, the shock front does not seem to have expanded much further than the central 111 kpc, or 15″, from the gas bridge center²⁰. This indicates that the shock front emerged rather recently at the gas bridge due to the recent collision of cluster material. Consequently, it is challenging to determine its direction of motion.

To characterize the speed of the shock front and estimate its Mach number ℳ, we use the Rankine–Hugoniot jump conditions for the density and temperature jumps (e.g., Landau & Lifshitz 1959; Russell et al. 2010), which respectively result in

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{dens}}={\left[\frac{2n_{\mathrm{e},1}/n_{\mathrm{e},2}}{\gamma +1-(\gamma -1)n_{\mathrm{e},1}/n_{\mathrm{e},2}}\right]}^{1/2}\end{array}$$(11)

and

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{temp}}={\left[\frac{{(\gamma +1)}^2(T_{\mathrm{X},1}/T_{\mathrm{X},2}-1)}{2\gamma (\gamma -1)}\right]}^{1/2}.\end{array}$$(12)

In the above equations, n_e, 1/T_X, 1 and n_e, 2/T_X, 2 are the gas density/temperature of the central box bin (where the shock front is located) and the pre-shock area respectively. Moreover, γ = 5/3 is the adiabatic index for a monatomic gas. As the pre-shock area, we consider the regions covered by the two neighboring box bins to the bin where the n_e and T_X jumps occur. We find ℳ_dens = 1.09 ± 0.10 and ℳ_temp = 1.31 ± 0.43, with the two values being consistent within the uncertainties and indicating a shock front velocity u_sh, close to the speed of sound c_s = 962 ± 52 km s⁻¹. Specifically, for ℳ_dens and ℳ_temp we respectively obtain u_sh = 1050 ± 111 km s⁻¹ and u_sh = 1259 ± 419 km s⁻¹. The Mach number indicates a weak shock front. Assuming the two Mach numbers are independent, their weighted mean is ℳ_comb = 1.11 ± 0.10, consintent with being subsonic at ≈1σ. The 3σ upper limit for the Mach number is ℳ_comb ≤ 1.41.

Given that the shock front seems to be restricted in the central bin, we assume it has traveled ≲111 kpc since it was created. Assuming that the shock’s velocity is ∼1100 km s⁻¹, the age of the shock is ≲100 Myr. Furthermore, one can determine the kinetic energy E_kin that flows through the shock front as (e.g., Gu et al. 2019; Zhang et al. 2019; Omiya et al. 2023)

$$\begin{array}{c}\hfill E_{\mathrm{kin}}=\frac{1.17}{2}{m}_{\mathrm{p}}{n}_{\mathrm{e},\mathrm{pre}}{u}_{\mathrm{sh}}^{3}S(1-\frac{n_{\mathrm{e},1}^2}{n_{\mathrm{e},2}^2}),\end{array}$$(13)

where m_p is the proton mass and S is the 2D area of the shock, which we take to be S = 15″ × 90″ = 7.41 × 10⁴ kpc². Both the u_sh and S estimates have considerable uncertainties; thus the final E_kin estimate is mostly indicative and not a precise constraint. Nevertheless, we find that E_kin ∼ 2.4 × 10⁴⁴ ergs/s; this is a slightly below-average amount of energy for a supersonic shock front and is expected given the low Mach number.

Moreover, the velocity of the shock front is representative of the cluster collision velocity on the plane of the sky (e.g., Omiya et al. 2023), which seems to be of order ∼1000 − 1300 km s⁻¹. As explained before, the emergence of the shock also implies that the two clusters are almost at the same cosmic distance, that is, separated by < (R_500, west + R_500, east)≈1.4 Mpc. Such a difference in cosmic distance would correspond to a cosmological Δz_Hubble ≈ 4 × 10⁻⁴. Thus, we can safely assume that ≳87% of the observed Δz = 0.003 is due to their relative peculiar motion and their collision velocity in the line of sight is ∼520 km s⁻¹. The transverse and line of sight approximate collision velocities imply a total collision velocity of 1100 − 1400 km s⁻¹, with the axis of motion forming a ∼22° −27° angle with the plane of the sky.

6. Summary

Pre-merging pairs of galaxy clusters are of special interest, as they can provide valuable insights into the early phases of cluster interaction, the physical processes taking place during the colliding of the hot ionized gas, and the formed gas bridges as well as into the impact of the initial merging on the behavior of the individual clusters. In this work, we performed an in-depth X-ray analysis of SPT-CLJ2228-5828, a distant (z ≈ 0.77) pair of clusters that have started colliding. From past SZ and optical analyses, SPT-CLJ2228-5828 was considered to be a candidate for a dissociative post-merger where the gas was separated from the two distinctive galaxy overdensities and WL mass peaks. Our work reveals that this system is a pre-merging double cluster system instead, with the radii of the two clusters overlapping and the gas of their outskirts interacting. Our work comprises of two parts: the study of the undisturbed sides of the two galaxy clusters and the analysis of their interaction region and its gas bridge.

When one focuses on the clusters’ half-circle regions opposite the interaction axis, the two separate clusters show highly similar X-ray properties with T_X, CE ≈ 3.2 − 3.4 keV, L_X, CE ≈ 0.7 × 10⁴⁴ erg s⁻¹, M_gas ≈ (1.8 − 2.1)×10¹³ M_⊙, and Y_X, CE ≈ (5.7 − 7.0)×10¹³ M_⊙ keV. Using the Y_X mass proxy, we estimated the total cluster mass to be M_500, YX = 2.10 ± 0.20 × 10¹⁴ M_⊙ and M₅₀₀ = 2.41 ± 0.22 × 10¹⁴ M_⊙ for the eastern and western cluster, respectively. We also estimated the WL masses of both clusters. While the western cluster is well detected with an estimated M_500, WL, which is consistent with the western cluster mass determined by X-rays, the eastern cluster is only marginally detected by its WL signal. The total mass of the system as determined by X-rays and WL is respectively ∼38% and ∼30% higher than the past SPT-SZ results. Moreover, we also extract the SB, n_e, T_X, P, and K profiles of both clusters, which show a typical behavior. Interestingly, if one ignores the interacting side of the clusters due to the ongoing merging, both systems consistently follow well-established cluster scaling relations. They also appear to be relaxed, according to their X-ray concentration and central n_e values. All of these findings suggest that interacting clusters may be used as single objects for cosmological and astrophysical studies if one focuses on the cluster sides opposite the interaction region. Furthermore, by simulating the expected SPT-SZ image of the merger based on its observed X-ray profiles, we find that the true SPT-SZ image is consistent with the nature of the merger; thus, the two clusters are detected as one only due to the limited angular resolution of SPT. Finally, assuming the two clusters have the exact same cosmic distance, we conclude that they are colliding with ∼1100 − 1400 km s⁻¹ with a ∼22° −27° angle compared to the plane of the sky.

The region between the two merging clusters shows characteristics of great interest. Firstly, a significant gas overdensity in addition to the cluster outskirts is robustly detected, indicating the existence of a separate gas bridge between the two clusters with ≈333 kpc length. We computed its gas mass to be M_gas = (1.37 ± 0.34)×10¹² M_⊙. The baryon overdensity of the bridge (excluding the expected gas from the clusters) is δ_b ≈ 320 at its center and δ_b ∼ 75 − 100 at its edges. Overall, the studied bridge properties suggest that there was excess gas between the clusters before the collision started, for example, in the form of a short filament. A weak shock front found at the center of this gas bridge shows distinctive gas density and temperature jumps, up to n_e ≈ 10⁻³ cm⁻³ and T_X ≈ 5.5 keV, respectively. The Mach number is estimated to be ℳ ≈ 1.1. The shock front also exhibits a mild pressure and entropy increase compared to the expected value if no shock was present. The age of the shock is estimated to be ≲100 Myr, and its velocity is u_sh ∼ 1100 km s⁻¹. The kinetic energy that flows through the shock was found to be E_kin ∼ 2.4 × 10⁴⁴ ergs/s. This is one of the rare detections of shock fronts in interacting pre-merging clusters. Finally, SPT-CLJ2228-5828 is the only merging system that has been analyzed using X-ray data showing a gas bridge and shock front and is located at such a high-z.

Acknowledgments

We thank the anonymous referee for their helpful comments. This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with GO program 16488. This work is based on observations obtained at the international Gemini Observatory, a program of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation on behalf of the Gemini Observatory partnership: the U.S. National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This research makes use of data products derived from observations taken by the 3D-HST Treasury Program (HST-GO-12177 and HST-GO-12328) with the NASA/ESA Hubble Space Telescope. This study was also based on observations obtained by XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA member states and the USA (NASA). These observations were acquired through a successful GO proposal (PI: H. Zohren). K.M. acknowledges support in the form of the X-ray Oort Fellowship at Leiden Observatory. The Bonn and Innsbruck authors acknowledge support from German Federal Ministry for Economic Affairs and Energy (BMWi) provided through DLR under projects 50OR2002, 50OR2106, and 50OR2302, from the German Research Foundation (DFG) under grant 415537506, and the Austrian Research Promotion Agency (FFG) and the Federal Ministry of the Republic of Austria for Climate Action, Environment, Mobility, Innovation, and Technology (BMK) via grants 899537, 900565, and 911971. Support for Program number GO-16488 was provided through a grant from the STScI under NASA contract NAS5-26555.

References

All Tables

All Figures

Fig. 1. Fully corrected, background-subtracted XMM-Newton smoothed EPIC count-rate image of the SPT-CLJ2228-5828 cluster pre-merger in the 0.5 − 2 keV band. The R500 and R200 radii of each cluster are displayed with white dashed and cyan circles respectively. The region between the solid green circles represents the area of the sky background estimation. The emission limit of the nearby active galactic nucleus and its host galaxy is also marked with a white dashed circle.

In the text

Fig. 2. Overlay of the signal-to-noise ratio map of the WL mass reconstruction (shown as white contours starting at level S/N = 1.5 in steps of ΔS/N = 0.5) and a color image of the inner cluster regions, employing the HST F140W, F814W, and F606W images for the R, G, B channels. For the non-parametric mass reconstruction shown in this image, the positions of the WL centers were left free to vary. In contrast, for the mass constraints presented in Sect. 3.2.4 we fix them to the X-ray cluster centers. The two positions, as well as the BCG positions, agree well within ≲1σ. The cyan circle indicates the SZ centre from Bleem et al. (2015), while the magenta star marks the BCG of the western cluster component. The red contours show the smoothed X-ray MOS count-rate in steps of 1.385 × 10−5 cts.

In the text

Fig. 3. Spectroscopic redshift information of the two clusters. Left: Spatial distribution of the Gemini-studied galaxies with spectroscopic z. The color bar represents the z value of each galaxy. The light green and purple circles represent the 1′ area around the eastern and western clusters (X-ray peaks), respectively, and display the regions within which galaxies are assumed to belong to the respective cluster. Galaxies with bright yellow or deep blue colors do not belong to the cluster system. Middle: Spectroscopic galaxy z histogram for the eastern cluster. The black line represents the X-ray z normalized likelihood. Right: Same as in middle panel but for the western cluster.

In the text

Fig. 4. Surface brightness profiles (PIB-subtracted) for the eastern and western clusters (left and right panels, respectively) in the 0.5 − 2 keV band. The total (cluster+CXB) PSF-convolved profiles are displayed in black, while the cluster-only PSF-deconvolved profiles are shown in the filled blue band, which represents the 68.3% error band. The R500 and R200 values are displayed with the dashed green and blue vertical lines respectively. The CXB level is displayed with the black horizontal line.

In the text

Fig. 5. X-ray radial profiles for the eastern and western clusters (green and blue colors, respectively). The R500 values are displayed with the dashed vertical lines. Upper left: Gas (electron number) density ne profiles with their best-fit models. The R200 values are also displayed in this panel. Upper right: Temperature profiles (projected, 2D) with their best-fit deprojected, 3D models. Due to the limited data, no model uncertainties could be derived. Bottom left: Inferred pressure profile best-fit models computed as P = ne × TX, 3D. The 1σ error band is underestimated since it accounts only for the uncertainty of the ne profile assuming no uncertainty of the TX, 3D profile. Bottom right: Inferred entropy profile best-fit models computed as $K=T_{\mathrm{X},3\mathrm{D}}\times n_{\mathrm{e}}^{-\frac{2}{3}}$. The 1σ error band is underestimated since it accounts only for the uncertainty of the ne profile assuming no uncertainty of the TX, 3D profile.

In the text

	Fig. 6. Image from SPTpol of the cluster merger (top) and the expected SPT-SZ Compton-y image based on the derived X-ray ne and TX profiles and the SPT-SZ angular resolution (bottom). The two images are consistent. That is, the SZ data are consistent with the X-ray analysis of this work.
In the text

	Fig. 7. Zoomed-in real (top) and simulated (bottom) image of the system in the 0.5 − 2 keV band. The simulated image assumes only the emission of the clusters and no gas bridge. The displayed boxes represent the regions in which the SB and ne in Figs. 8 and 9 were measured.
In the text

Fig. 8. Surface brightness profile (CXB-subtracted) across the merging axis for the boxes shown in Fig. 7, measured in the 0.5 − 2 keV band. The black points represent the measured SB values. The gray ones represent the expected SB values from cluster projection effects and assuming no excess gas bridge. The position of the eastern and western X-ray peaks are displayed by the dashed green and blue vertical lines respectively. The horizontal black dashed line shows the CXB level. The SB ratio between the measured values and the expected values if there was no gas bridge is shown in the bottom panel.

In the text

Fig. 9. X-ray profiles across the cluster merging axis for the temperature (top left), gas (electron number) density (top right), pressure (bottom left), and entropy (bottom right). The TX, P, and K profiles are measured for larger bins than the ones in Fig. 7, with 31.5″ width as explained in Sect. 5.2. The gray points in the ne profile describe the expected ne value due to projection effects if no gas bridge was present. All X-ray properties are elevated compared to the case where no gas bridge or shock front exists between the two clusters.

In the text