Observations of a Pre-Merger Shock in Colliding Clusters of Galaxies

TL;DR

This paper reports the discovery of a unique pre-merger shock in a galaxy cluster pair, providing new insights into early cluster collision dynamics and the thermal history of large-scale structures.

Contribution

It presents the first observation of a pre-merger shock propagating outward along the equatorial plane, an early-phase phenomenon not previously documented.

Findings

Discovery of a shock in a pre-merger galaxy cluster pair.

Shock propagates outward along the equatorial plane, unlike typical merger shocks.

Supports models predicting shock energy dissipation outside clusters.

Abstract

Clusters of galaxies are the largest known gravitationally-bound structures in the Universe. When clusters collide, they create merger shocks on cosmological scales, which transform most of the kinetic energy carried by the cluster gaseous halos into heat. Observations of merger shocks provide key information of the merger dynamics, and enable insights into the formation and thermal history of the large-scale structures. Nearly all of the merger shocks are found in systems where the clusters have already collided, knowledge of shocks in the pre-merger phase is a crucial missing ingredient. Here we report on the discovery of a unique shock in a cluster pair 1E 2216 and 1E 2215. The two clusters are observed at an early phase of major merger. Contrary to all the known merger shocks observed ubiquitously on merger axes, the new shock propagates outward along the equatorial plane of the merger. This discovery uncovers an important epoch in the formation of massive clusters, when the rapid approach of the cluster pair leads to strong compression of gas along the merger axis. Current theoretical models predict that the bulk of the shock energy might be dissipated outside the clusters, and eventually turn into heat of the pristine gas in the circum-cluster space.

Tables

Table 1. Table 1: List of Chandra , XMM-Newton , and Suzaku observations used in this paper.

Instrument	Observation	Date	Initial exposure (ks)	Clean exposure (ks)
Chandra ACIS	20778	2018-07-23	38.0	37.6
Chandra ACIS	21131	2018-07-25	32.1	31.7
Chandra ACIS	21132	2018-07-26	36.1	35.6
Chandra ACIS	21133	2018-07-27	25.6	25.3
Chandra ACIS	21134	2018-07-29	17.1	16.9
XMM-Newton EPIC	0800380101	2017-11-11	139.0	79.5 (MOS), 53.7 (pn)
Suzaku XIS	807085010	2012-05-16	28.6	24.9
Suzaku XIS	807084010	2012-11-29	15.8	11.2

Equations

$$T(x,y)=\sum _{\mathrm{k}}{G}_{\mathrm{k}}(R)T(x_{\mathrm{k}},y_{\mathrm{k}})/\sum _{\mathrm{k}}{G}_{\mathrm{k}}(R),$$

$$S\approx \Lambda (T,A)\int n_{\mathrm{e}}^{2}dl,$$

$$P=T{S}^{1/2},$$

$$K=T{S}^{-1/3}.$$

$$\mathcal{M}=\sqrt[]{\frac{(\gamma +1)\frac{P_{\mathrm{post}}}{P_{\mathrm{pre}}}+(\gamma -1)}{2\gamma }}\approx \mathrm{1.5.}$$

$$n(R)=\left\{\begin{array}[]{c l}Cn_{0}R^{-\alpha_{1}}&R\leq R_{\rm edge},\\ n_{0}R^{-\alpha_{2}}&R>R_{\rm edge},\end{array}\right.$$

$$\mathcal{M}=\sqrt[]{\frac{2C}{\gamma +1-(\gamma -1)C}}=\sqrt[]{\frac{3C}{4-C}},$$

$$\delta \mathcal{M}=\frac{6}{(4-C{)}^2\mathcal{M}}\delta C$$

$$\frac{T_{\mathrm{post}}}{T_{\mathrm{pre}}}=\frac{[2\gamma {\mathcal{M}}^2-(\gamma -1)][(\gamma -1){\mathcal{M}}^2+2]}{(\gamma +1{)}^2{\mathcal{M}}^2}=\frac{5{\mathcal{M}}^4+14{\mathcal{M}}^2-3}{16{\mathcal{M}}^2},$$

$$\mathcal{M}=\sqrt[]{\frac{\left(8\frac{T_{\mathrm{post}}}{T_{\mathrm{pre}}}-7\right)+\sqrt[]{{\left(8\frac{T_{\mathrm{post}}}{T_{\mathrm{pre}}}-7\right)}^2+15}}{5}},$$

$$\delta \mathcal{M}=\frac{4\mathcal{M}}{\sqrt[]{{\left(8\frac{T_{\mathrm{post}}}{T_{\mathrm{pre}}}-7\right)}^2+15}}\delta \left(\frac{T_{\mathrm{post}}}{T_{\mathrm{pre}}}\right).$$

$$n_{\mathrm{e},\mathrm{post}}^{1-\gamma }{T}_{\mathrm{post}}=n_{\mathrm{e},\mathrm{pre}}^{1-\gamma }{T}_{\mathrm{pre}},$$

$$t_{\mathrm{eq}}\approx 2\times 10^9{\left(\frac{T}{10^8\mathrm{K}}\right)}^{3/2}{\left(\frac{n_{\mathrm{p}}}{10^{-4}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{\mathrm{ln\Lambda }}{40}\right)}^{-1}\mathrm{yr},$$

$$\Delta F_{\mathrm{K}}=\frac{1}{2}{\rho }_{\mathrm{pre}}{v}_{\mathrm{s}}^3\left(1-\frac{1}{C^2}\right),$$

Full text

• 1

RIKEN High Energy Astrophysics Laboratory, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

• 2

SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, the Netherlands

• 3

ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA, Dwingeloo, The Netherlands

• 4

Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, The Netherlands

• 5

International Centre for Radio Astronomy Research, Curtin University, GPO Box U1987, Perth, WA 6845, Australia

• 6

Hamburger Sternwarte, Universität Hamburg, Go jenbergsweg 112, 21029, Hamburg, Germany

• 7

MTA-Eötvös University Lendület Hot Universe Research Group, Pázmány Péter sétány 1/A, Budapest, 1117, Hungary

• 8

Institute of Physics, Eötvös University, Pázmány Péter sétány 1/A, Budapest, 1117, Hungary

• 9

Department of Physics, University of Strathclyde, Glasgow, G4 0NG, UK

• 10

Raman Research Institute, C. V. Raman Avenue, Sadashivanagar, Bangalore 560080, India

• 11

Department of Physics & Electronics, Rhodes University, PO Box 94, Grahamstown, 6140, South Africa

Here we present a study of an pre-merger cluster pair 1E 2216.0-0401 and 1E 2215.7-0404 at a redshift of 0.09, primarily with 149/139/44 kiloseconds of Chandra/XMM-Newton/Suzaku X-ray observations. As shown in Figure 1, the cores of the two clusters are separated by a distance of 640 kpc. The central regions of the two clusters have similar ICM temperatures ($4.7\pm 0.1$ keV and $5.4\pm 0.1$ keV), suggesting a nearly equal-mass merger. Most interestingly, the Chandra image reveals a wedge-like feature extended towards the southeast on the equatorial plane between the two clusters (Figure 2). The wedge is clearly seen in the X-ray residual map where the X-ray haloes of the two clusters are modeled (Supplementary Figure 1) and subtracted, suggesting that it cannot just be the overlapping of the two clusters. The leading edge of the wedge is located about 450 kpc away from the collision axis. A radial profile of X-ray surface brightness as a function of distance from the collision axis shows a clear discontinuity across the edge (Supplementary Figures 2-4), indicating a sharp drop in gas density. The surface brightness steepening at the discontinuity significantly (with 98% confidence) deviates from a smooth King profile or an exponential profile.

The equatorial plane is found to be hotter than the surrounding area in the projected temperature map made with the XMM-Newton data (Figure 3 and Supplementary Figure 6). The peak temperature, $\sim 8.1\pm 1.4$ keV, is found at the leading edge of the wedge. A secondary peak is seen at the north side of the equatorial plane. No hot structure is seen at the other directions around the two clusters, suggesting that the observed temperature increase cannot be merely a temperature gradient internal to each of the clusters. Combining the temperature and surface brightness distribution around the wedge, we find a significant pressure discontinuity at its leading edge. This strongly favours the presence of a shock, rather than a cold front, at the location of the observed surface brightness discontinuity. An X-ray entropy map reveals the signs of heating extending towards a larger distance than the observed shock front, hinting for previous equatorial shocks at an earlier stage of the merger (Supplementary Figure 7).

We use different approaches to calculate the velocity of the observed shock. By measuring the surface brightness discontinuity, the ICM behind the shock front is found to be compressed by a mean value of $2.1\pm 0.5$, giving a shock Mach number $\mathcal{M}=1.8\pm 0.5$. A better way to estimate the shock speed is to determine the temperature jump at the shock (Supplementary Figure 5), because this method is less affected by the unknown shock geometry [17, 18, 19, 20]. The average temperature jump obtained from the different instruments weighted by their errors is $1.6\pm 0.4$, giving a shock Mach number $\mathcal{M}=1.6\pm 0.3$. The observed temperature jump is nearly consistent with the value expected from adiabatic compression, suggesting that the ICM heating is mostly done with the shock compression. Assuming that the shock is propagating from the collision axis with an increasing velocity [16], the shock has an age $>250$ Myr. Comparing it to the age of the possible axial shock ($\sim 50-100$ Myr) reported in [18], the equatorial shock seems to have formed substantially earlier.

Numerical simulations [15, 16] predict that merger shocks along the equatorial plane should appear in pairs, however in our target only one is visible in the X-ray image (Figure 2). We then searched for signatures of shocks using 144 MHz radio observations from the LOw Frequency ARray (LOFAR), and 235 MHz and 325 MHz observations from the Giant Metrewave Radio Telescope (GMRT). The radio images reveal a complex radio structure at the northwest of the equatorial plane, where a second equatorial shock is expected. The radio feature is mainly dominated by a pair of radio tails of $\sim 330$ kpc size (source A in Figure 2), and a bright patch with a narrow extension of $\sim 450$ kpc towards northwest (source B). A close inspection of the LOFAR and GMRT images shows that the radio tails are probably associated with a reported cluster member galaxy (Supplementary Figure 8; [21]). The 610 MHz GMRT image reveals a compact radio core in the probable host galaxy. Therefore, we conclude that source A is associated with a member AGN, but the nature of source B still remains unclear.

The spectral measurement between 144 MHz and 325 MHz reveals a steep spectrum of the radio tails (source A), with an average spectral index of $-2.3\pm 0.1$. Naturally, the radio spectrum is expected to steepen away from the AGN due to ageing effect, however, the observed spectral index distribution (Figure 3 and Supplementary Figure 9) appears to be too flat to be explained by radio radiation. The cooling cosmic ray particles must have been re-energized, probably through a Fermi-type shock acceleration process [22, 23]. The re-acceleration balances out the radiative loss, inducing a flattening of the spectral index at the shock [11]. Furthermore, source A seems to overlap, in projection, with the heated ICM structure ($k$T $\geq 7.0$ keV) at the northwest of the equatorial plane (Figure 3). It is therefore likely that the re-acceleration of the radio tails and the ICM heating are related, probably powered by the same shock (a twin to the southeast shock) or an adiabatic compression triggered by the merger.

Contrary to the north region, we do not detect significant radio emission at the southeast shock. If the north radio sources are revived from cooling cosmic ray particles, we may conclude that the bulk of such particles is not present at the southeast. The direct injection of relativistic particles out of the thermal pool might not be efficient, as the observed shock speed, $\mathcal{M}\sim 1.6$, is smaller than minimal Mach number required for acceleration in several theoretical models [24, 25, 26].

The X-ray observations of the sharp surface brightness discontinuity and the associated high temperature structure provide the best evidence to date for the presence of an equatorial shock in an early stage of on-axis merger. The powerful shock at southeast contains a kinetic energy flow rate of $\sim 8.1\times 10^{44}$ erg s*-1*, about two times of the combined X-ray luminosity of the two clusters. The northwest shock, if exists, carries a similar amount of kinetic energy as the southeast one.

Our observation resembles the results from the state-of-the-art theoretical models. The latest hydrodynamic codes [15, 16] predicted that the very first pair of merger shocks should occur before the core collision, and propagate outward along the equatorial plane. These shocks are formed by the interaction of the outer regions of the two clusters. Previously, there is only tentative observational evidence for this theory, including the X-ray temperature increase between the early-stage mergers A85 [27], A115 [19], A399-A401 [28], and A1758 [29], as well as a Sunyaev-Zel’dolvich (SZ) decrement reported by [30] between the two sub-clusters of CIZA J2242+5301. However, A85 and A115 are off-axis mergers, the high temperature regions could be explained alternatively by a side-by-side interaction of regular shocks driven by the infalling subclusters. For the on-axis mergers A399-A401 and A1758, no sharp X-ray surface brightness discontinuity could be found with the available data. Hence, it is unclear whether the observed temperature breaks are caused by shocks, particularly equatorial shocks, or by a milder adiabatic compression of the two merger clumps during the early merger stage. Similarly, adiabatic compression could also explain the SZ feature in CIZA J2242+5301. Moreover, CIZA J2242+5301 is already at a late stage of merger. It is hard to detect a clean feature from its early stage as it would have already been mixed. We believe that the detection of the equatorial shock in 1E 2216.0-0401 and 1E 2215.7-0404 provides the ideal and most direct support for the theoretical model.

The discovery of equatorial shock completes the inventory of merger shocks. It reveals an important merging epoch right before the core collision, when shocks are formed at the contact interface and propagate vertically to the merger axis. For the first time, we measure the energy release by the equatorial shock, which may prove to be one of the important dissipation processes of major mergers. As shown in the numerical simulations [15, 16], the equatorial shocks have the highest velocity among merger shocks, allowing them to reach large distance ($\geq 3$ Mpc) before they disappear. These shocks, together with some axial shocks that form later, might dump a significant portion of kinetic energy beyond the cluster boundary, providing a feasible way to heat the inter/circum-galactic medium in the cluster surroundings.

Methods

We assume that $\rm H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{\rm M}=0.27$ and $\Omega_{\rm\Lambda}=0.73$. At the distance of the 1E2216 cluster, 1*′* corresponds to 91.8 kpc. All spectra are fitted using the C-statistics, and the quoted uncertainties are 1$\sigma$.

Chandra data processing

Chandra observed 1E2216 on 2018 July $23-29$ for a total exposure of 149 ks with the ACIS-I detector, pointing at the bridge region between the two merging clusters. The datasets are reprocessed using the software chandra_repro in the latest CIAO 4.10 package. Light curves are extracted using 200 s time bins in the $0.5-12.0$ keV for the point-source-free region. These light curves show no flares, therefore the entire observing time period is used. For each observation the bright point sources are removed using the tool celldetect with the detecting threshold of 3$\sigma$ on the $0.5-8.0$ keV image. The image is further checked by eye for fainter point sources.

The X-ray background is removed using the blank sky background dataset. For each observation, the tool blanksky is run to locate the blank sky files, and to reproject the blank sky to match the observation. The background files are reprocessed to have the same gain, timing, and bad pixel properties as the observation data.

The exposure map is created for each observation using the tool mkexpmap in the $0.5-8.0$ keV band. To take into account the energy dependence, we provide a model spectrum to mkexpmap, which has a temperature of 5 keV, abundance of 0.3 solar, and a Galactic absorption column of $7.3\times 10^{20}$ cm*-2* [1]. The edges of the detector are excluded from the following analysis, as they show some possible artefacts caused by the exposure correction. The final flux image is then created by removing the background from the original image, and dividing it by the exposure map. Figure 2 shows the flux image in the $0.5-8.0$ keV. It is smoothed by a Gaussian kernel with an adaptive scale to include about 50 net counts in each kernel.

To highlight the possible surface brightness excess on the equatorial plane of the merger, we subtract the X-ray emission from the main haloes of the two clusters. The classical beta profile is used to model the projected 2-D gas distribution in the $0.5-8.0$ keV band for each cluster. The equatorial regions, as shown in Supplementary Figure 2(a), are excluded in the fits. The shape of the King profile, the central position and the ellipticity of the X-ray haloes, are all set free in the fits. The best-fit model is then smoothed by a Gaussian kernel with the same scale as is done to the flux image. As shown in Supplementary Figure 1, the resulting beta-like image is used to calculate the residual map (Figure 2).

For the spectral analysis, the source spectral files are extracted separately from each of the observations, and the responses and ancillary response matrices are created using mkacisrmf and mkwarf, respectively. Background spectra are extracted from the blank sky datasets. The spectra of the same region are fit simultaneously between 0.5 keV and 7.0 keV using the SPEX package with SPEXACT 3.04. The raw spectra are re-grouped with the optimal binning method provided by the SPEX package.

To obtain an overview of the two clusters, we analyse two spectra, each covering the central 350 kpc region of one cluster. The Chandra spectra are fitted with a single redshifted cie model for the emission measure, temperature, average metallicity, and the redshift. The proto-Solar abundances from [2] are adopted. A foreground Galactic absorption column of $N_{\rm H}=7.3\times 10^{20}$ cm*-2* is included in the fits. The best-fit temperatures are $4.6\pm 0.1$ keV and $5.4\pm 0.1$ keV, the average metallicities are $0.36\pm 0.06$ solar and $0.28\pm 0.05$ solar, and the redshifts are $0.095\pm 0.004$ and $0.091\pm 0.004$, for the northeast and southwest clusters, respectively. The obtained parameters agree well with those reported in the previous Suzaku work. Based on the mass-temperature relation found in [3], the approximate total masses of the two clusters are $2.5\times 10^{14}$ solar mass and $3.2\times 10^{14}$ solar mass.

XMM-Newton data processing

XMM-Newton observed 1E2216 on 2017 November 11-13 for a total exposure of 139 ks. The latest XMM science analysis system (SAS) and the extended source analysis software (ESAS) are used to process and calibrate the data obtained with the European Photon Imaging Camera (EPIC) onboard XMM. The tools emchain and epchain are used to create the reduced EPIC-MOS and pn event files, respectively. Light curves are extracted in 100 s bins and screened for background flares with the mos-filter and pn-filter tasks. The final exposures reduce to $\sim 80$ ks for the MOS and $\sim 54$ ks for the pn data. The events with PATTERN values greater than 12 (MOS) / 4 (pn) and non-zero FLAG values are also excluded. The exposure-corrected images are made in the 0.5 to 2 keV band using a bin size of 40 pixels. The MOS image is used to identify point sources from the diffuse components. The sources detected above a flux threshold of $5\times 10^{-14}$ ergs cm*-2* s*-1* are excluded in the analysis. The Chandra image is used to cross-check the point source detection made with the MOS.

XMM-Newton provides a better statistics of the target cluster than Chandra and Suzaku. The source spectra, particle background, and response files are prepared by the mos-spectra and pn-spectra tools. The remaining sky background, including the foreground Galactic soft X-ray emission and the unresolved cosmic X-ray background (CXB) are taken into account as spectral components in the fits. The Galactic emission is modelled by two thermal plasma components (the SPEX cie model) with fixed temperatures of 0.08 keV and 0.3 keV. The 0.08 keV component is free of absorption while the 0.3 keV one is absorbed by the Galactic ISM. The abundances and redshifts of the two components are fixed to unity and zero, respectively, while their normalizations are set free to vary. The CXB component is represented as an absorbed power law, with a fixed photon index of 1.4 and flux of $6.4\times 10^{-8}$ ergs cm*-2* s*-1* sr*-1* [4].

We examine the XMM-Newton spectra from the ICM properties in the central 350 kpc regions of the two clusters. The spectra are re-grouped with the optimal binning method. The ICM component is modelled by one redshifted cie component, absorbed by the Galactic neutral materials. The column density is fixed to $7.3\times 10^{20}$ cm*-2*. The temperature, average metallicity, redshift, and emission measure of the thermal component are left free in the fits. The best-fit temperatures are $4.7\pm 0.1$ keV and $5.1\pm 0.1$ keV, and the average metallicities are $0.42\pm 0.04$ solar and $0.33\pm 0.03$ solar, for the northeast and southwest clusters, respectively. The redshifts of the two objects are $0.094\pm 0.001$ and $0.090\pm 0.001$. These values are consistent with those determined with the Chandra spectra.

Suzaku data processing

The Suzaku X-ray Imaging Spectrometer (XIS) observed 1E2216.0-0401 and 1E2215.7-0404 on May and November 2012 for exposures of 28.6 and 15.8 ks, respectively. Data reduction is performed with HEAsoft version 6.24 and the latest CALDB version (20160607). To suppress instrumental background, an event screening with the cosmic-ray cutoff rigidity COR2 $>$ 6 GV is applied. Further reductions are applied to the XIS1 to mitigate background increase due to the change of the charge injection (http://www.astro.isas.jaxa.jp/suzaku/analysis/xis/xis1_ci_6_nxb/). Because of the damage from a micro-meteoroid strike, XIS2 is not in use. The resultant clean exposures are about 25 and 11 ks for 1E2216.0-0401 and 1E2215.7-0404, respectively.

For the spectral analysis, we model the sky background in the same way as is done for XMM-Newton. All the XIS spectra are fit simultaneously.

Radio data preparation

Radio observations at 235, 325, 610 and 1300 MHz were obtained with the GMRT under project codes 28_087 and 31_079. At 235 and 610 MHz (dual-frequency observing), the target was observed on 2015-08-19 for 2.2 hours over an effective bandwidth of 18 and 33 MHz, respectively. At 325 MHz, the target was observed on 2016-12-24 for 6.1 hours over an effective bandwidth of 33 MHz. At 1300 MHz, the target was observed on 2016-12-23 for 2.8 hours over an effective bandwidth of 33 MHz. All data was processed using the SPAM pipeline, which included flux scale and bandpass calibrations using 3C 48, and multiple rounds of direction-dependent calibration and imaging. The final images have sensitivities of 570, 130, 57, and 96 $\mu$Jy/beam at 235, 325, 610, and 1300 MHz, respectively. The resolutions are 14.8” x 9.6”, 11.7” x 8.0”, 6.5” x 3.8”, and 2.4” x 1.9”, respectively. The 1300 MHz image is not used in this paper, because it is not deep enough to detect the main part of the target source.

The LOFAR high band antenna (HBA) data were obtained as part of projects LC7_003 and LC9_014 and the target was observed for 4 hours on four separate occasions between 2017-01-19 and 2017-11-21 with a frequency coverage of 120-168 MHz. All the data were passed through the standard direction independent calibration pipeline, and made use of the LOFAR Long Term Archive (LTA) compute facilities for efficient processing. The data were then processed through the latest version of the LOFAR Two-metre Sky Survey (LoTSS) direction dependent calibration pipeline [5]. This pipeline uses kms [6] to calibrate for ionospheric effects and errors in the beam model and ddfacet [7] to apply the derived solutions during the imaging. To further enhance the quality of the images a post processing step was added to remove all sources from the field besides the target and to self calibrate the data to improve the accuracy of the calibration solutions at that location. The final data has a sensitivity of 0.27 mJy/beam in the vicinity of the target when imaged with a resolution of 17.6” x 6.3”.

Thermal structure of the pre-merger system

To obtain a global view of the thermal structure, we map the two-dimensional distribution of the ICM temperature in the early merger cluster. The map shown in Figure 3 is made based on the method described in [8]. First, we randomly select a set of pixels to be the centers of extraction regions. The distribution of “knots” is sufficiently dense, with a separation $<10^{\prime\prime}$ between any two adjacent knots, to ensure that the map is fully covered by the extraction regions. For each knot, a spectrum is extracted from a circular region containing a net counts of about 7000 (or a signal to noise of 84) in the $0.5-7.0$ keV band. The background is corrected with the XMM-Newton blanksky. Typical radii of the extraction regions are $\sim 0.3^{\prime}-1.5^{\prime}$. The best-fit temperatures, determined from the fits with the SPEX cie model, are assigned to the corresponding knots. The temperatures of the remaining pixels are obtained through interpolation. For a pixel at position $(x,y)$, we define a circular region containing 7000 counts, with a radius of $r(x,y)$. The temperature at $(x,y)$ is then calculated through a weighted average over all the knots $(x_{\rm k},y_{\rm k})$ in the circular region,

[TABLE]

where $R$ is the distance from $(x,y)$ to $(x_{\rm k},y_{\rm k})$, and $G_{\rm k}$ is the Gaussian kernel with the scale parameter $\sigma$ fixed to $r(x_{\rm k},y_{\rm k})$. The use of Gaussian function produces a smooth temperature map, whose resolution is set by the extraction region. As the main body of the map is essentially derived through interpolation, it is little affected by the bias from the possible poor fits in individual regions.

The statistical errors on the temperature are calculated with the error command in the fits. In the fits of a few regions, the error calculation could lead to a smaller C-statistic than the original fits, we update the fits with the new value and run the error calculation again. The temperature errors map is derived through the same technique as the temperature map. As shown in Figure 3, the temperature errors are typically $\sim 0.3$ keV for the cluster central regions, $\sim 0.6$ keV for the bridge region between the two clusters, and $\sim 0.8$ keV at the southeast shock region. The errors are about 10% of the values in the temperature map, except for the outermost regions at $r\geq 7^{\prime}$, where the errors exceed 20% of the local temperatures.

As a crosscheck, we further calculate the temperature maps using the standard weighted Voronoi Tesselations (WVT) algorithm [9, 10, 11, 12, 13] (http://www.phy.ohio.edu/~diehl/WVT/) and the adaptive circular binning (ACB) [14, 15]. The WVT binning divides the map into 150 non-overlapping regions, each contains 7000 clean counts in $0.5-7.0$ keV. The temperatures between two adjacent WVT bins are determined independently. The ACB method extracts spectra from a circular region, which contains 7000 clean counts, centered on each “pixel” as shown in Supplementary Figure 6(c). The ACB regions significantly overlap with each other in the cluster outskirt. As can be seen in Supplementary Figure 6, the three temperature maps obtained through different methods agree well with each other. For most of the regions, the discrepancies among the maps are found to be $\leq 0.5$ keV, smaller than the typical statistical errors from the spectral fits.

As shown in Figure 3, the peak of ICM temperature can be seen at the southeast wedge, probably due to the heating by the equatorial shock. The temperature decreases, from $\sim 8$ keV to $\sim 6$ keV, across the surface brightness discontinuity. It seems that the temperature gradient in the map is slightly milder than the deprojected temperature profile for the post-shock and pre-shock regions, which reveals a jump from $8.1$ keV to $5.2$ keV. This is probably because that hot spot in the post-shock region has been blurred by the interpolation method (Eq. 1 in the Method section).

Combining the Chandra image and the XMM-Newton temperature map, we can further obtain the approximate distributions of gas pressure and entropy. The electron number density distribution can be directly determined from the X-ray surface brightness $S$, as

[TABLE]

where $\Lambda(T,A)$ is the cooling function that depends on the ICM temperature and metallicity, and the electron density $n_{\rm e}$ is integrated along the line of sight. Ignoring the cooling function will introduce a uncertainty of about 10% in the pressure and entropy maps, which is similar to the typical error of the temperature map.

By approximating electron density as $S^{1/2}$, the pseudo-pressure and pseudo-entropy are then defined as

[TABLE]

[TABLE]

The Chandra $0.5-2.0$ keV image is used for the pressure and entropy maps. The Chandra image is smoothed with the XMM-Newton point spread function, and then convolved with the the same adaptive smoothing kernel $r(x,y)$ as the XMM-Newton temperature map. A region of both high pressure and entropy might indicate the presence of a shock. As shown in Figure 3 and Supplementary Figure 7, the southeast shock is indeed shown as a high pressure and entropy region, while such a dramatic structure is not visible at the northwest. The high pressure structure nicely traces the southeast shock, and shows an abrupt drop outside the front. Reading directly from the map, the ratio of the post-shock pressure $P_{\rm post}$ to pre-shock pressure $P_{\rm pre}$ is $\sim 2.4$, and the shock Mach number can be estimated as

[TABLE]

It agrees with the values obtained from the surface brightness and temperature jump analysis. Such a gas pressure jump can potentially be confirmed with a dedicated high spatial resolution observation of the Sunyaev-Zel’dovich effect [16, 17, 18, 19]. The pressure map also reveals a mild gradient at the north, extending up to $\sim 3^{\prime}$ from the collision axis.

As shown in Supplementary Figure 7, the high entropy feature at the southeast appears to extend to a larger radius than the detected shock front, even taking into account the relatively low resolution of the temperature map. This might indicate the presence of other ICM heating processes in the past, such as previous equatorial shocks and/or adiabatic compression by the merger. The high entropy is also seen along the equatorial direction at the north of the collision axis, although the average gas entropy in the north is lower than the value in the south.

Data availability

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.