© The Authors 2025

1. Introduction

Current cosmological models and numerical simulations predict that the majority of the missing baryons in our Universe reside in faint galaxy cluster outskirts and interconnecting filaments of the cosmic web (Cen & Ostriker 1999; Davé et al. 2001). This warm-hot intergalactic medium (WHIM) contributes to the growth of galaxy clusters through a process of slow, constant accretion, as the infalling filamentary gas virialises within their gravitational potential wells.

Only a small number of X-ray bridges and filaments have been studied so far. This is partially because the diffuse WHIM is still challenging to detect with today’s instruments. Numerical cosmological simulations predict WHIM temperatures of 10⁵ to 10⁷ K and densities of 10⁻⁷ to 10⁻⁴ cm⁻³ (Cen & Ostriker 2006; Haider et al. 2016). In other words, features in the WHIM are inherently characterised by soft X-ray emission, low surface brightness, and poor signal to noise ratios (S/Ns) due to dominant contributions to the signal from both astrophysical and instrumental backgrounds. It is clear that a great deal of care needs to be taken when modelling the spectra of these complex regions to account for all of the different contributions from other emission components.

Current observational information on filaments is derived from either absorption or emission. Observations in emission have so far been limited to a small number of systems, MACS J0717.5+3745 (Ebeling et al. 2004), Abell 399/401 (Sakelliou & Ponman 2004), Abell 222/223 (Werner et al. 2008), Abell 2811 (+offset)/2804/2801 (Sato et al. 2010), Abell 3558/3556 (Mitsuishi et al. 2012), Abell 2744 (Eckert et al. 2015), Abell 3391/3395 (Alvarez et al. 2018; Reiprich et al. 2021; Veronica et al. 2022, 2024), Abell 2029/2033 (Mirakhor et al. 2022), Abell 98N/98S (Sarkar et al. 2022), and Abell 3667/3651 (Dietl et al. 2024), as well as A3530/32 and A3528-N/S (Migkas et al. 2025). Temperature measurements were obtained only for five of these systems. For two of the five, it has been reported that the filament emission is dominated by the emission of the intracluster medium (ICM), as demonstrated by the excessively large temperatures. There is an ongoing debate about the origin of these X-ray filaments. While some of the detected bridges could indeed be classified as WHIM filaments, others are due to other physical processes, such as shock-related compression heating between the two merging atmospheres or emission of ram-pressure stripped tails from infalling halos.

One of the most extensively studied merging galaxy clusters in almost every available wavelength is MACS J0717.5+3745 (RA 07h17m32.1, Dec +37°45′21″), an extremely massive, intermediate-redshift (z = 0.5458) Hubble Frontier Fields cluster (Ebeling et al. 2001). MACS J0717.5+3745 is also one of the most complex galaxy clusters known to date, which has hosted at least four different sub-cluster collisions. The collective mass from the many ingested dark matter halos makes this a great target for gravitational lensing studies and additional exploration of local substructures via the Sunyaev-Zeldovich (SZ) effect due to the high temperatures. The mergers have also resulted in a complex ICM morphology, characterised by discontinuities, such as shock fronts and cold fronts, as well as many peculiar features in the X-ray and radio bands. In addition to the gravitational lensing study, Jauzac et al. (2018) also looked at many of these X-ray substructures in detail. One of its most intriguing features is the prominent X-ray bridge located in the S-SE of the cluster, which was previously studied by van Weeren et al. (2016, 2017) using Chandra data.

This paper aims to investigate the X-ray bridge in detail, statistically disentangling the diverse distribution of complex emission signatures associated with the filament using Bayesian X-ray analysis (BXA) methods and nested models, along with a new dynamically-adaptive method for the full instrumental background modelling of each region for both Chandra and XMM-Newton. This paper also presents new thermodynamical maps using a joint modelling method that combines the available data from Chandra, XMM-Newton, and ROSAT, as well as trend-modeled residual maps where the average cluster thermodynamical properties over radial distances have been removed. This paper additionally presents previously unpublished data from 90 GHz MUSTANG2 observations of the SZ effect in the direction of MACS J0717.5+3745. In Section 2, we explain the methods used in the data reduction, followed by Section 3, which explains the instrumental and X-ray background model and analysis of both Chandra and XMM-Newton. Section 4 describes the main results of the thermodynamic maps and the properties of the filament, while Section 5 discusses these results. Throughout the paper, we assume a standard Λ cold dark matter cosmology with Ω_m = 0.286, Ω_Λ = 0.714, and H₀ = 69.6 (Bennett et al. 2014). Consequently, at z = 0.5458, 1 arcmin corresponds to 387.42 kpc. The wilms abundance table is adopted for all plasma emission and photoelectric absorption models in the discussed spectral models (Wilms et al. 2000). Unless stated otherwise, the error bars correspond to a 68% confidence interval.

2. Observations and data reduction

2.1. Chandra observations

The data from four Chandra pointings (see Table 1) were reprocessed using repro from the level-1 event lists with the standard software packages using the most recent versions of CIAO and CALDB (versions 4.16.0 and 4.11.0). The good time intervals (GTIs) of the observations after filtering periods of flaring are summarised in Table 1. All the observations were used for creating the surface brightness images; however, only observations 4200, 16235, and 16305 were considered in the further spectral analysis, as observation 1655 experienced residual soft proton flaring contamination which would bias spectral analysis. Combining all the available Chandra observations gives approximately 228.9 ks of observing time; however, with the removal of observation 1655, only 211.8 ks were used for spectral analysis. Because of the low S/N, blank-sky background files were explicitly not used in favour of using epoch-specific stowed background files normalised by the data over background count ratio in the 10-12 keV energy range. When dealing with extremely faint X-ray substructures, the stowed background files give additional control over the systematics in the background. Merged broadband images were produced in the 0.5 to 7 keV band, shown in Fig. 1(left). the spectra of all regions were extracted from each Chandra observation using the specextract routine in the latest version of CIAO (version 4.16.0).

Fig. 1. Images of the full FOV of Chandra (left) and XMM-Newton (right) of MACS J0717.5+3745. The surface brightness excess to the west/south-west are several foreground structures in between us and the cluster, indicated by the ‘foreground emission’ label. The XMM-Newton image is over-saturated in the cluster centre so that this can be seen more clearly. The green dashed circle indicates the approximate R200 of the cluster. Intensity scale is in flux per pixel (0.492 × 0.492 arcsec for Chandra, and 4.1 × 4.1 arcsec for XMM pn).

2.2. XMM-Newton observations

Complementary to the Chandra data, we also retrieved archival XMM-Newton observations of MACS J0717.5+3745, as listed in Table 1. Three pointings showing sufficient data quality were downloaded and reduced using both the XMM Science Analysis System (SAS) software (v20.0.0) and the current calibration files (CCFs) by following the procedure detailed in the Extended Source Analysis Software (ESAS v20.0.0) cookbook. To get cleaned event files for each detector, we used the standard emchain and mos-filter routines for the metal oxide semi-conductor (MOS) data, and the epchain and pn-filter routines for the pn data. mos-filter and pn-filter remove periods of anomalously high count rate by calculating GTIs through the espfilt routine. The total GTI for each observation are also shown in Table 1.

Images were extracted for each detector of each observation using the ESAS tasks mos-spectra and pn-spectra, while background images were created using the ESAS tasks mos_back and pn_back, before being combined using the comb routine. Point sources were detected using a wavelet detection algorithm, wavdetect, by first making a psfmap for XMM-Newton with a constant size of 9 arcsec. This is based on the estimate that at 1.5 keV, the 1-σ integrated volume of a 2D Gaussian for each pixel ranges from 7 arcseconds on-axis to around 11 arcseconds at the edge of the FOV and our selected box region of interest is near the central pointing. Merged broadband images from all observations were produced in the 0.5 to 7 keV band using the merge_comp_xmm routine, shown in Fig. 1(right). A zoomed-in, over-saturated view of the cluster centre and filamentary structure (with pointsources removed) can be seen in Fig. 2.

Fig. 2. Chandra (left) and XMM-Newton (right) images of MACS J0717.5+3745, restricted to the same FOV as the SZ data and thermodynamical maps, visible as the box in Fig. 1. Intensity scale is in flux per pixel (0.492 × 0.492 arcsec for Chandra, and 4.1 × 4.1 arcsec for XMM pn). Part of the XMM-Newton image is over-saturated so that the filament structure is more visible. The white lines indicate the regions used for modelling the spectra of the outskirts of the cluster and the group together with the filament. The green dashed circle indicates the approximate R200 of the cluster.

The spectra of all regions are extracted from each XMM-Newton observation using the standard ESAS tools mos-spectra, mos_back, pn-spectra and pn_back. The XMM-Newton regions are first converted into observation and instrument-specific detector coordinates; however, because of the many point sources and complex region geometries, the complete region expression becomes too long for the fits data subspace. Therefore it is necessary to first turn each sky region minus any overlapping point sources into a region file in detector coordinates for extraction with ESAS.

All spectra from both Chandra and XMM-Newton were optimally binned using the method described in Kaastra & Bleeker (2016) and all astrophysical parameters were linked through all spectra for a joint fit. Unless otherwise specified, all spectra in this work were produced in a similar way. These regions were finally fit using SHERPA version 4.16.0 (from CIAO-4.16), using the XSPEC version 12.13.1e model library and AtomDB version 3.0.9.

2.3. CLASH weak-lensing observations

In this study, we used ground-based weak-lensing data products from the CLASH program, as presented in Umetsu et al. (2014). These products were derived from deep BVR_Ci′z′ imaging with Subaru/Suprime-Cam, complemented by UV imaging from Megaprime/MegaCam and near-infrared (NIR) observations from WIRCAM on the Canada–France–Hawaii Telescope (see also Medezinski et al. 2013). For further details on data reduction and weak-lensing analysis, we refer to Umetsu et al. (2014).

2.4. MUSTANG2 Sunyaev-Zeldovich observations

MACS J0717.5+3745 was observed by the 100-meter Green Bank Telescope (GBT) under project codes AGBT17A_340 and AGBT17B_266 (in semesters 2017A & 2017B, respectively) using the MUSTANG-2 instrument, a 215-element continuum bolometer array operating at 90 GHz (Dicker et al. 2014). The observations were conducted using Lissajous daisy scans with radii of 3′ and 3$\stackrel{\prime }{.}$5 (Romero et al. 2020, 2023). The data processing procedure was identical to that presented by Romero et al. (2020), amongst others. The resulting map, produced using the MUSTANG-2 pipeline, is shown in the lower left panel of Fig. 11. The total integration time on target was 2.6 hours and the resulting map has a resolution of 10″, with an RMS noise level of ∼28 μJy/bm. Besides the astrophysical contaminants and the instrumental noise, MUSTANG2 data is affected by large-scale filtering of the astrophysical signal, which affects the imaging quality.

3. Background analysis

3.1. Chandra stowed instrument background

The Chandra stowed background files are generated in a similar procedure as the blank-sky background files and they are scaled by the data over background count ratio in the 10–12 keV energy range. The Chandra instrumental background can be described as a mixture of various signal contributions across different energy bands. Specifically, cosmic rays hit the detectors from every direction and create a continuum emission. Additionally, the particle interactions can produce instrumental X-ray fluorescence lines from the various elements in and around the detector.

Of the available background datasets, the total stowed background merged exposure time was 2.91 Ms for the year 2000, 2.57 Ms for 2005, and 1.68 Ms for 2009. The spectra were extracted from the full field of view (FOV) of each available year and were modelled following a similar process as described in Siemiginowska et al. (2010) for Chandra and Mernier et al. (2015) for XMM-Newton, with independent fits for each of the detectors. Special care was also taken regarding the faint and very faint imaging modes, resulting in modified models that involve these different cases. The stowed background data was modelled in the 0.7 keV to 10 keV band for both the ACIS-I and ACIS-S detectors. The continuum emission is defined by a broken power law for the faint imaging submode and by a power law in the very faint imaging submode, while all instrumental lines were modelled using Gaussians and frozen to their best fits. Bartalucci et al. (2014) noted that the shape of the particle background is fairly consistent, and only the normalizations in different epochs are changing from year to year; therefore, we refrained from including redundant information for all years. We simply give detailed results of the continuum emission and instrumental lines in Table 2. Elements corresponding to line energies used in our background model were documented in other works (Bearden & Burr 1967; Suzuki et al. 2021).

The line energies that did not exactly correspond to reported lines are marked in the table with ‘’, while the 2.6 keV line is marked with ‘*’ since it was not previously mentioned in other reference materials. We note that this line showed up only in the ACIS-I stowed backgrounds and is unlikely to be an instrumental line; rather, it must be some artifact of the response itself. We also note that it was included in the ACIS-I model as an empirical improvement for completeness. All Chandra spectral fits reported in this paper were performed in the 0.7 to 7.0 keV range.

3.2. XMM-Newton instrumental background

Similarly to the Chandra stowed backgrounds, the XMM-Newton filter wheel closed instrumental backgrounds can also be described as a mixture of various spectral contributions across different energy bands. Besides the ‘hard particle’ (HP) or quiescent particle background (QPB) components shared with Chandra, XMM-Newton cluster source observations can additionally suffer from a residual contribution from a ‘soft-particle’ (SP) component; these are remnants of the quiescent soft proton contamination, which may still pervade throughout the data despite having been excised from the observation during the GTI cleaning step. These ‘soft-protons’ are particles trapped in Earth’s magnetosphere, which are accumulated and concentrated by the flight optics during the orbit cycle, affecting around 40% of all observation time.

The hard particle components were modelled using all currently available FWC data, corresponding to a total merged exposure time of around 4.77 Ms. Specifically, around 1.958 Ms of data from MOS1, 1.924 Ms of data from MOS2, and 0.8841 Ms of Extended Full-Frame pn data. Spectra were extracted from the full FOV of all available FWC observations, and modelled following a similar process as described in Mernier et al. (2015), also with independent fits for each of the detectors.

The FWC data were modelled in the 0.3 keV to 10 keV band for the MOS1 and MOS2 detectors and in the 0.4 keV to 10 keV range for pn. The continuum emission is fitted by a broken power law, while all elemental lines were modelled using Gaussians described in Table 3. Because these HP components are noise contributions directly from the detectors and are unrelated to the mirrors, these instrumental background models are not convolved by the ancillary response file (ARF).

The residual SP components can be also described as an ARF-unfolded power-law component with an index limited to a range between 0.1 and 1.4, which affects source spectra. However, the inclusion of the additional power law did not contribute to an improvement in the background model fit statistics, indicating that the soft protons were fairly well treated in the data reduction.

To further check for any residual soft proton contamination, we compared the area-corrected count rates between the ’inside’ and ‘outside’ regions of the FOV for each detector in each observation, shown in Table 4 (De Luca & Molendi 2004). All of the diagnostic values are under the lowest soft proton threshold of 1.15, indicating that none of the event lists are additionally contaminated by soft protons. As such, in this paper, the soft proton parameters were removed from the background model, which reduced the overall complexity. All XMM-Newton spectral fits reported in this paper were performed in the 0.45 to 7.0 keV range for MOS and 0.3 to 7.0 keV for pn.

3.3. The astrophysical X-ray background

With the instrumental backgrounds accounted for, special attention must be given to carefully model the astrophysical components of the background before we can accurately constrain any of the physical parameters of the cluster component of the signal. It is especially in the faint signal-to-noise ratio (S/N) and signal-to-background (S/b) regimes, such as cluster outskirts or filamentary structures, the parameters are all very sensitive, so any small changes to the instrumental background can affect the astrophysical background, which can both collectively bias the final results.

The cosmic X-ray background (CXB) includes the unresolved point sources, distant active galactic nuclei (AGNs), typically modelled with a photoelectrically absorbed power law with its slope frozen at 1.41 and a free normalization (Hickox & Markevitch 2006). The Galactic foreground sources can be described as a combination of several plasmas in collisional ionisation equilibrium (apec). They include the local hot bubble (LHB), the galactic halo (GH), and one more, recently discovered component, interpreted as a super-virialised (SV) component of our galactic halo, sometimes referred to as the ‘galactic corona’. The LHB is an approximately 0.1 keV plasma originating from a local supernova remnant that we are currently residing in. It can be effectively modelled using the ROSAT All Sky Survey (RASS) data with the support of the soft X-ray bands from Chandra and XMM-Newton. Additionally, our galactic halo and its supervirial component are described by photoelectrically absorbed (tbabs, wilms), thermal APEC plasma models (Wilms et al. 2000). For all of these local, foreground CXB components, the redshift was frozen to 0 and the abundances frozen to solar, with the temperatures and normalizations left free to fit. Unless otherwise stated, the 1.41 power-law photon index was kept frozen. Meanwhile, the neutral hydrogen column density for the CXB and in all subsequent fits was frozen to N_H = 8.36 × 10²⁰ cm⁻², corresponding to the sum of the atomic and molecular hydrogen column densities along the line of sight (LOS)¹ (Kalberla et al. 2005; Willingale et al. 2013).

To constrain the local CXB parameters, we used a carefully constructed ‘blank field’ for our observations by using all available data outside of the cluster R₂₀₀, from both the Chandra and XMM-Newton data, as well as spectra from the ROSAT All Sky Survey. The ROSAT All Sky Survey spectrum was extracted from an annulus with radii 0.15 and 1 degrees around the cluster centre, similar to van Weeren et al. (2017). We use a fit range for ROSAT PSPC-C/B between 0.09 and 2.0 keV based on previous calibration database plots of on-axis effective area curves. Meanwhile, the Chandra and XMM-Newton observations utilised the entire FOV of the observations with most point sources and emission regions removed. Additionally, the Chandra data excluded bright edges and the XMM-Newton data also had approximately 1 arcminute from the edges of the observations removed, as can be seen in Fig. 3.

Fig. 3. Emission-free regions taken from the observations to fit the CXB components, with Chandra on the left, and XMM-Newton on the right. The box is the FOV used for Fig. 2 and the thermodynamical maps, corresponding to the total FOV of the SZ data. The green dashed circle indicates the approximate R200 of the cluster. The excluded circles are removed sources, the large excluded box region to the west corresponds to removed foreground galaxies, and another smaller box removes the emission of the filamentary structure.

Using the available data from each instrument, the CXB was constrained using a total of 13 spectral data groups. The best-fit parameters of the foreground components are often degenerate with each other. For this reason the LHB, GH, SV, and unresolved point sources were determined using generous priors with the Bayesian X-ray analysis tool, which uses a nested sampling algorithm (UltraNest²) for a Bayesian parameter estimation and model comparison employed as a tool for faint S/N data within the AGN community (Buchner 2021; Buchner et al. 2014).

The top part of Table G.1 shows the results of the different CXB components, with and without the inclusion of the ROSAT All Sky Survey data, using both a background subtracted and background modelled fits. Without the ROSAT data, the LHB temperature was fixed to 0.09 keV and a free normalization. Using all available data and the background models, the best-fit results using all 13 spectra gave a LHB temperature of 0.108 ± 0.02 keV, GH temperature of 0.156 ± 0.002 keV, and SV temperature of 0.697 ± 0.011 keV, with a SV normalization around 10% that of the GH, consistent with other studies.

In all of the spectral analyses described in this paper, the spectral data were never combined; instead, all the data were fit simultaneously. For each of the spectral data groups, a 20% systematic error was applied in Sherpa/Xspec to account for any cross-calibration related discrepancies in measurements between the instruments. A more detailed discussion regarding the cross-calibration can be found in Appendix B.

Once the CXB parameters are constrained, it is possible to more accurately constrain the cluster emission and, consequently, to create thermodynamical maps and trend-averaged maps of MACS J0717.5+3745 to better study and visualise the ICM physics and internal structures. Specifically, we used a modified version of the contour binning algorithm from Sanders (2006), which we first translated into Python from C++, and have released as pycontbin³. This contour binning algorithm creates regions with a pre-selected S/N, while conserving the general contours of the data by grouping neighbouring pixels of similar surface brightness, and it yields statistically independent regions. An added benefit of this modified contour binning code is that it can additionally take a static value for a point spread function (PSF) as a constraint and merge map bins that are smaller than the given value, ensuring bins larger than the instrument PSF. The contour binning maps were derived from an exposure corrected XMM-Newton mosaic image, scaled to have the approximate number of counts per pixel as the sum of observations. The smallest bin in this map corresponds to a circle of around 10.6 arcseconds radius, around the size of the PSF of XMM-Newton. Like the CXB region, all contour bin map regions were converted into polygons and were extracted with both Chandra and XMM-Newton, with the idea that the Chandra resolution could be exploited in the joint fitting procedure as a correction to the fit results in the cases where the complex polygon region geometries occasionally fall under the size of the XMM-Newton psf. This resulted in a joint fit between 25 different spectra for each spatial bin, 13 from the CXB region, and 12 spectra from each spatial bin.

We again used the model for a plasma in collisional ionisation equilibrium (apec) with photoelectric absorption (tbabs), but now with fixed values for the metallicity at 0.3 Solar and redshift at 0.5458. We additionally scaled each region by its area in arcminute² using a model constant, so that we could simultaneously fit the large cxb region with each cluster spectra. The CXB parameters are first frozen to their best-fit values before fits are performed to constrain the cluster temperatures (kT) and normalizations (norm). Finally, the ICM metallicity is thawed to constrain its value for each region. This fitting procedure is repeated for every bin and for every instrument independently before being performed on the combined data from all observatories and again for the combined data from all available instruments. Unless otherwise mentioned, this spectral fitting methodology applies to any other spectral analysis considered for the rest of the paper.

4. Results

4.1. Thermodynamical maps

Exploiting this joint-fitting methodology using Chandra and XMM-Newton with the full instrumental background modelling resulted in around 110 spatial bins for the electron density map, around an order of magnitude improvement in the map resolution compared to the around 12 spatial bins shown in the previous map reported in van Weeren et al. (2017), which only exploited the available Chandra data and the blank-sky fields. Despite the smaller bin sizes, all map bins shown in this paper are larger than the XMM-Newton PSF, and our large reported temperature of ≥20 keV is consistent with the previous reported temperature in van Weeren et al. (2017).

The thermodynamical maps were computed from a scaled, exposure-corrected XMM-Newton mosaic image, with S/N values of 120, 100, 70, and 50 (14 400, 10 000, 4900 and 2500 counts per bin, respectively). Regions from S/N 120 were extracted from XMM-Newton, and the other bins from only Chandra to make use of the better spatial resolution. To compute asymmetric errors, we resample the multivariate distribution of the parameters from the covariance matrix provided by the fit results 100 000 times, using the median and the 16th and 84th percentiles of the resampled distribution for the 68% confidence intervals. Shown in Fig. 4 is the S/N 120 temperature map and S/N 100 density map, along with their combined S/N pressure and entropy maps, all of which had relative errors and propagated uncertainties lower than 20%. Combining the S/N 120 and 100 region maps yields a better-resolved gradient between map regions at the expense of having completely statistically independent bins. This is our favoured method for improving the map resolution without suffering from over-smoothing, as seen in other binning methods. We discuss the maps in more detail in Sect. 5.1.

Due to the nature of merging systems, projection effects may bias the results due to potential asymmetries along the projected LOS. Moreover, because this cluster is a highly-disturbed, non-spherically symmetric, merging cluster with possibly more than four distinct sub-cluster populations in different merging states, it is difficult to easily approximate the integration volume due to uncertainties in the projected geometry. To begin with a simplification, we first assume a flat, pancake-like approximation for the LOS depth of 1 Mpc for each map pixel in each bin, shown in Fig. 4. The units given in the thermodynamic maps in Fig. 4 are, cm⁻³×(l/1Mpc)^−1/2 for the density, keVcm⁻³×(l/1Mpc)^−1/2 for the pressure, and keVcm²× (l/1 Mpc)^1/3 for the entropy.

Fig. 4. Maps of temperature (top left), density (top right), pressure (bottom left), and entropy (bottom right) with corresponding units of keV, cm−3×(l/1Mpc)−1/2, keVcm−3×(l/1Mpc)−1/2, and keVcm2× (l/1 Mpc)1/3. The temperature and electron density maps were created using binning with a S/N of 120 (14400 counts) and 100 (10000 counts), respectively, while the pressure and entropy maps were produced using the results from the top two S/N maps. Density maps were derived from the normalization fits assuming a 1 Mpc projected LOS depth. Spectral fits were made using data from both Chandra and XMM-Newton for the temperature, and only Chandra for the density, each using the CXB best-fit results from ROSAT, Chandra, and XMM-Newton. Fits for these maps were performed using the most appropriate background model for each region according to the AIC.

The temperature, kT, and electron density, n_e, were derived directly from the spectral fitting, although the electron density was determined from a normalization parameter, norm, formally defined as

$$\begin{array}{c}\hfill \mathrm{norm}=\frac{{10}^{-14}}{4\pi {(d_a(z+1))}^2}{\displaystyle \int n_{\mathrm{e}}}{n}_{\mathrm{p}}\mathrm{d}V,\end{array}$$(1)

where d_A is the angular diameter distance, z is the redshift, and ∫n_en_p is the integrated cluster emissivity over the chosen volume of gas in the cluster. The electron density n_e and ion density n_p are related as n_e = 1.18n_p. The electron pressure and entropy are then determined as P_e = n_ekT and $K=kT{n}_{\mathrm{e}}^{-2/3}$, respectively.

Given the inherent unknowns and uncertainties associated with the true three-dimensional geometry of MACS J0717.5+3745, we explore the systematic effects resulting from the simplifying assumption of a constant 1 Mpc LOS depth in our original thermodynamic maps. To quantify potential biases introduced by this approximation, we implement a geometrical correction derived from a model of the cluster in a 3D box. Specifically, we generated an effective line-of-sight length map using five different 3D spherical beta models parameterised according to the lensing-detected sub-clusters listed in Table 7, with an additional beta model added to include the group at the end of the filament. The best-fit results are given in Table 5. The emissivity distribution was then integrated from the centre of the 3D simulated cube until 50% of the total flux along the LOS was reached, shown in Fig. 5.

Fig. 5. Map of the line-of-sight effective length derived from the 3D numerical simulation. The five sub-cluster centres are shown as white circles, corresponding to the different beta models shown in Table 5.

Furthermore, we quantified the effect of the geometric correction used in this paper with additional maps of the fractional change and associated sigma change, using the propagated uncertainties derived from spectral fitting for each region, which can be seen in Appendix E (Fig. E.1) and Appendix F (Fig. F.1).

4.2. Azimuthal trend-divided maps

We additionally created trend-divided maps to enhance the small azimuthal variations present in the cluster shown in Fig. 7. We began by calculating a centroid that is the average between the centre of the central contour bin (Bin 0), the position of the peak X-ray surface brightness pixel, and the best-fit results from a beta model fit to the cluster centre. We then constructed a scatter plot of the radial distribution for each physical quantity from the cluster centroid and fit for the average trend using a function expressed as

$$\begin{array}{c}\hfill f(r)=A{(1+{(r/B)}^2)}^{(-3C/2)}{(1+{(r/D)}^2)}^{(-3E/2)},\end{array}$$(2)

Fig. 6. Maps of temperature (top-left), density (top-right), pressure (bottom-left), and entropy (bottom-right) with corresponding units of keV, cm−3, keVcm−3, and keVcm2. The temperature and electron density maps were created using binning with a S/N of 120 (14 400 counts) and 100 (10 000 counts), respectively, while the pressure and entropy maps were produced using the results from the top two S/N maps. Density maps were derived from the normalization fits assuming an effective length determined from the best-fit results of five beta models, shown in Table 5. Spectral fits were made using data from both Chandra and XMM-Newton for the temperature, and only Chandra for the density, each using the CXB best-fit results from ROSAT, Chandra, and XMM-Newton. Fits for these maps were performed using the most appropriate background model for each region according to the AIC.

Fig. 7. Geometry-corrected, trend-divided maps of temperature (top-left), density (top-right), pressure (bottom-left), and entropy (bottom-right).

where r is the radial distance from the centroid, and A, B, C, D, and E are free parameters, following the method described in Ichinohe et al. (2015). Since we wanted to remove the general thermodynamical trends from the cluster, we did not use the full annuli. Instead, we restricted the radial distribution of physical quantities to the angles corresponding to the same region used to fit for the cluster outskirts. Trend-divided maps are then created for each physical quantity by calculating a residual for every bin, between the original value and the model-predicted value at the given radius. Since these trend-divided maps are unit-less, they are likewise less sensitive to geometric uncertainties imposed by our chosen cluster geometry.

4.3. Spectroscopic disentanglement of the filament signal

The spectroscopic fitting of the filamentary region was performed in several steps. We kept the standard assumption that it is a faint signal contaminated by instrumental backgrounds, as well as both astrophysical foreground and background signals; however, we extended the astrophysical description of the background to also include the local atmospheres of nearby systems. The filament lies between MACS J0717.5+3745 and a group to the southeast, inside their respective R₂₀₀ emission areas; thus, to properly constrain the temperature and density of the filament, it is first necessary to disentangle the different contributions to the complex signal.

To constrain the cluster and group-related contributions to the filament signal, large arc regions were extracted, approximately at the same radial distance at which the filament structure is displaced from the cluster and group centres, inside, but close to R₂₀₀. These regions were then independently fit using a similar method described in Sect. 4.1, using CXB parameter constraints as determined from the optimal CXB fit and then applying these posteriors as priors to the subsequent fits. The fit results of the CXB + Cluster and CXB + Group regions are shown in the centre of Appendix G (see Table G.1).

As we did not know the exact multi-phase composition of the filamentary gas, we tried to statistically disentangle the multi-phase nature of the region using Bayesian X-ray analysis and model selection using the Akaike and Bayes information criteria (AIC and BIC) described in Appendix A. Table 6 shows the AIC and BIC values when comparing a combination of different spectral models that could effectively describe the filament region. To fit all of the various models, we used only the CXB and filament spectra to fit the temperature and normalization of the filament data, while all other parameters were frozen to their best fits. In these data, there were 3116 data points, and approximately 3114 degrees of freedom, so the AIC and BIC values are fairly close to each other.

We first compared a simple CXB + Filament (apec) model with a CXB + Filament (gadem) model, which applies a Gaussian distribution for the emission measure versus temperature to an apec model. Despite being often used to fit multi-phase plasmas, gadem was punished by the statistical measure for its extra complexity compared to apec and it was rejected as an option for the more complex model comparisons. We then proceeded to fit a CXB + cluster + filament model and a CXB + group + filament model, before finally fitting a rather complex model with CXB + cluster + group + filament. This complex model involved the joint-fitting of 49 independent spectra: 13 CXB-associated spectra and 12 spectra for each of the cluster, group, and filament regions shown in Fig. 2.

Based on the different criteria, all of the various apec models were close enough in their likelihoods with each other that they were not outright rejected and they could theoretically could each be effectively used to describe the data. That being said, the most complicated of the models, which contained the CXB + cluster + group + filament, ended up with the best-fit results, despite being punished for its extra complexity. Interesting to note is that the ‘second-best’ model according to the AIC is the CXB + cluster + filament, and not the expected CXB + group + filament model.

Figures 8, 9, and 10 show the source and background spectral data from Chandra, XMM-Newton MOS1 and MOS2, and XMM-Newton pn of the filament region. Besides the full astrophysical and instrumental background models, each plot also shows the individual contributions to the total model from each astrophysical temperature component, along with a residual plot for the spectral data and the full model, and a second residual plot for the instrumental background model. The counts for each instrument are averaged across each energy bin and then binned for improved visualization. As can be seen in the Figures, the contribution from each temperature component was successfully constrained by the fitting. The XMM-Newton instrumental background lines between 1.2 and 1.75 keV were not well fit with either a Gaussian, Lorentzian, or Voigt profile, so we ignored that spectral energy range in the spectral fitting. Other features present in the residual were due to the selected background model (in the case of bknpow with no Gaussians) or the cross-calibration between the XMM-Newton pn, XMM-Newton MOS, and Chandra data.

Fig. 8. Chandra filament spectra, including all astrophysical and instrumental model components. Spectral data was averaged for each energy bin in all relevant observations, and then re-binned for better visibility. The middle and lower plots correspond, respectively, to the residuals for the astrophysical source data with respect to the full model and the instrumental background data and related instrumental background model.

Fig. 9. XMM-Newton MOS1+MOS2 filament spectra, including all model components. Spectral data was averaged for each energy bin in all relevant observations and then re-binned for better visibility. The middle and lower plots correspond, respectively, to the residuals for the astrophysical source data with respect to the full model, as well as the instrumental background data and related instrumental background model.

Fig. 10. XMM-Newton pn filament spectra, including all model components. The spectral data was averaged for each energy bin in all relevant observations, and then re-binned for better visibility. The middle and lower plots correspond, respectively, to the residuals for the astrophysical source data with respect to the full model and the instrumental background data and related instrumental background model.

4.4. The Sunyaev-Zeldovich effect

Traditionally, the ICM has been studied by using X-ray observations; on the other hand, millimeter-wave astronomy provides an exceptional observational tool to directly probe gas pressure via the SZ effect first reported by Sunyaev & Zeldovich (1972). The SZ effect is a spectral distortion in the frequency spectrum of the cosmic microwave background (CMB) radiation due to inverse Compton scattering by a thermal plasma. (for a recent review, see Mroczkowski et al. 2019). In practice, it is commonly observed as either a decrement (or increment) in the CMB intensity at frequencies under (or over) 218 GHz along lines of sight through clusters of galaxies, with a linear dependence on the electron density (as opposed to X-ray observations which are dependent on electron density squared). The change in the effective CMB temperature is proportional to the Compton-y parameter, which depends on the Thomson scattering optical depth, τ_e, and temperature of the hot electron gas, T_e, expressed as

$$\begin{array}{c}\hfill y_{\mathit{tSZ}}\equiv {\displaystyle \int \frac{k_{\mathrm{B}}{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2}}\mathrm{d}{\tau }_{\mathrm{e}}={\displaystyle \int \frac{k_{\mathrm{B}}{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2}}{n}_{\mathrm{e}}{\sigma }_T\mathrm{d}l=\frac{{\sigma }_T}{m_{\mathrm{e}}{c}^2}{\displaystyle \int P_{\mathrm{e}}\mathrm{d}l,}\end{array}$$(3)

where the Thompson cross section, P_e = n_ek_BT_e, is the pressure due to the electrons integrated over the LOS distance (Mroczkowski et al. 2019).

There can additionally be another spectral distortion to the CMB spectrum due to the Doppler effect, known as the kinetic SZ effect, caused by the cluster bulk velocity along the LOS on the scattered CMB photons and described by

$$\begin{array}{c}\hfill y_{\mathit{kSZ}}={\sigma }_T{\displaystyle \int \frac{-v_z}{c}}{n}_e\mathrm{d}l\equiv \frac{-v_{\mathrm{z}}}{c}{\tau }_{\mathrm{e}}.\end{array}$$(4)

where v_z is the cluster bulk velocity (positive for receding clusters), and τ_e is the electron optical depth. The observed change of specific intensity, ΔI_ν, can be expressed with respect to the CMB intensity I₀ as a function of the observed frequency, ν, given by

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{I}_{\nu }}{I_0}=f(\nu ,T_e)y_{\mathit{tSZ}}+g(\nu ,T_e,v_z)t_{\mathit{kSZ}},\end{array}$$(5)

where $I_0=\frac{2{(k_B{T}_{\mathrm{CMB}})}^3}{{(hc)}^2}=270.33\left[\frac{T_{\mathrm{CMB}}}{2.7255\mathrm{K}}\right]\mathrm{MJy}/\mathrm{sr}$, and f(x, T_e) describes the characteristic tSZ effect frequency dependence, and g(ν, T_e, v_z) describes the spectral dependence of the kSZ effect in the non-relativistic regime, given by Birkinshaw (1999) as

$$\begin{array}{c}\hfill f(\nu ,T_e)=-\frac{x^4{e}^x}{{(e^x-1)}^2}(xcoth\left(\frac{x}{2}\right)-4)(1+{\delta }_{\mathrm{tSZ}}(x,T_e)),\end{array}$$(6)

$$\begin{array}{c}\hfill g(x,v_z,T_e)=\frac{x^4{e}^x}{{(e^x-1)}^2}(1+{\delta }_{\mathrm{kSZ}}(x,v_z,T_e)),\end{array}$$(7)

where $x=\frac{h\nu }{k_B{T}_{\mathrm{CMB}}}\approx \nu /56.8\mathrm{GHz}$ is the dimensionless frequency, h is the Planck constant, ν is the observed frequency, k_B is the Boltzmann constant, and T_CMB = 2.7255K is the CMB temperature. The δ_tSZ(x, T_e) and δ_kSZ(x, v_z, T_e) terms correspond to the relativistic corrections for the tSZ and kSZ effects, and are dependent on the observing frequency and T_e, and in the kSZ case, also the LOS velocity. The relativistic corrections for the tSZ are computed using the work of Itoh & Nozawa (2003), while the kSZ corrections are computed using the analytical formula in Nozawa et al. (2006). Practically, these relativistic corrections are computed from tabulated values using SZpack, whereby the corrections are on the 10% level for these 18 keV temperatures (Chluba et al. 2012, 2013).

Here, we introduce the previously unpublished MUSTANG2 images in Fig. 11 and Fig. C. The MUSTANG2 data was corrected for the kinetic and relativistic SZ effect following a constructed model based on the best-fit results of the ‘F2’ model from Adam et al. (2017b); this model construction is described in detail in Sect. 4.5.

Fig. 11. Comparison plots between the geometry-corrected X-ray pressure map (top left), the geometry-corrected trend-divided X-ray pressure map (top-right), followed by SZ MUSTANG2 data (bottom-left), next to a kinetic and relativistic SZ corrected version of the data, shown here in KelvinCMB units (bottom-right). The white contours are from the radio images seen in van Weeren et al. (2017) and the green contours corresponds to the MUSTANG2 data in the bottom-left. The circles represent the strong lensing sub-cluster centres given in Limousin et al. (2016), which were used to reproduce the kSZ+rSZ model used in Adam et al. (2017b), displayed in Table 7. The full FOV of the MUSTANG2 data can be seen in Appendix C.

4.5. Correction of kinetic SZ effects via ICM physical modeling

Mroczkowski et al. (2012) and Sayers et al. (2013) both explored MACS J0717.5+3745 in the context of the SZ effect, making it the first discovery of a cluster system with a kinetic SZ component. Adam et al. (2017a) later produced spatially resolved kinetic SZ maps of MACS J0717.5+3745 with the New IRAM KIDs Array (NIKA) on the IRAM 30m telescope, followed by a joint X-ray/SZ study of the projected gas temperature distribution (Adam et al. 2017b). Recently, Adam et al. (2025) also used NIKA data to probe the turbulent gas motions in MACS J0717.5+3745, finding a kinetic to kinetic plus thermal pressure fraction P_kin/P_kin + th of around 20% and an energy injection scale of around 800 kpc.

Because MUSTANG2 is dominated by the tSZ and limited to a single frequency ∼90GHz, it does not have much constraining power for the kSZ signal. As we can not constrain our own kSZ model, we instead opt to build a kSZ model for MUSTANG2. To make the correction, we use the same model parameters as ‘F2’ constrained in Adam et al. (2017b), which used the NIKA data and extra information from X-ray imaging to break degeneracies between v_z and τ_e. We start with the fact that previous strong lensing observations in Limousin et al. (2016) have shown that there are at least four sub-clusters present in MACS J0717.5+3745, listed in Table 7, and shown as circles in Fig. 11. Due to the linear dependence of n_e in the SZ data, the 3D gas distribution is built as a sum of each independent sub-cluster, where each sub-cluster is modeled with a spherically symmetric β-model,

$$\begin{array}{c}\hfill n_e(r)=n_{e0}{[1+{\left(\frac{r}{r_c}\right)}^2]}^{-3\beta /2},\end{array}$$(8)

and n_e0 is the central electron density, r_c is the core radius, and β is the slope of the radial profile (Cavaliere & Fusco-Femiano 1978). Each sub-cluster additionally is assumed to have isothermal gas temperatures, determined using 30 arcsecond regions from XMM-Newton in Adam et al. (2017b), and constant bulk velocities, all given in Table 7. This gas distribution is then analytically integrated along the LOS for each sub-cluster as

$$\begin{array}{c}\hfill {\tau }_e\equiv {\sigma }_T{\displaystyle {\int }_{-\infty }^{+\infty }}{n}_{e}dl=\sqrt{\pi }{\sigma }_T{n}_{e0}{r}_c\frac{\mathrm{\Gamma }(\frac{3}{2}\beta -\frac{1}{2})}{\mathrm{\Gamma }\left(\frac{3}{2}\beta \right)}{[1+{\left(\frac{R}{r_c}\right)}^2]}^{\frac{1}{2}-\frac{3}{2}\beta },\end{array}$$(9)

where Γ represents the Gamma function, and R is the projected radius from the cluster centre. This can then be used to create the tSZ and kSZ models used to compare with the MUSTANG2 surface brightness maps, given by

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }{I}_{\nu }^{\mathrm{model}}}{I_0}& ={\sigma }_T{\displaystyle \sum _i}{f}_{\nu }(T_x^{(i)})\frac{k_B{T}_x^{(i)}}{m_e{c}^2}{\displaystyle \int n_e^{(i)}}\mathrm{d}l\hfill \\ \hfill & +{\sigma }_T{\displaystyle \sum _i}{g}_{\nu }(T_x^{(i)})\frac{v_z^{(i)}}{c}{\displaystyle \int n_e^{(i)}}\mathrm{d}l,\hfill \end{array}$$(10)

where T_x and v_z are scalar quantities for each subcluster.

Fig. 11 shows a comparison plot between a zoomed in view of the geometry-corrected X-ray pressure and the geometry-corrected, trend-divided pressure maps (top), followed by SZ MUSTANG2 data and a kSZ and rSZ corrected version of the data. The white box corresponds to the FOV of the MUSTANG2 data and it is the same as the thermodynamical maps shown in Figs. 6 and 7, while the new smaller green box shows the zoomed-in region in the two top plots. The white contours are the radio images in van Weeren et al. (2017), and the circles represent the sub-cluster centres given in Table 5 that were used for the geometric correction of the maps and those given in Table 7 that were used in the kSZ model.

There is a relatively close agreement between the X-ray pressure peak and the SZ peak, with an offset corresponding to around 0.25 arcminutes, or around 100 kpc at this redshift, but this estimate is biased due to the resolution of the MUSTANG2 data and the size and shape of the corresponding bin from the X-ray pressure map. As can be seen from the trend-divided pressure map, there is also a clear offset between the radio relic emission and the X-ray pressure peak, likely due to the nature of the integrated X-ray emission and projection effects mentioned in van Weeren et al. (2017). Despite the large uncertainties and differences in the assumptions associated with the kSZ model, our results are consistent with the previously published results in Adam et al. (2017b). Considering the noise in the MUSTANG2 data, shown in Appendix C, the kSZ correction applies the most significant correction around sub-cluster B of around 2.9σ. The correction is less significant at the other sub-clusters, specifically –0.31σ, –0.96σ, and –0.28σ for sub-clusters A, C, and D, respectively.

As mentioned previously, the present work does not explore pressure fluctuations (see Khatri & Gaspari 2016; Romero et al. 2023, 2024), but we note that deeper SZ observations with, for example, MUSTANG2, could improve upon the constraints for MACS J0717.5+3745 presented in Adam et al. (2025).

5. Discussion

5.1. Thermodynamical maps

The temperature map on the top left of Fig. 6, shows an exceptionally hot cluster centre, at around 24 ± 4 keV, with a temperature dropping to 9.8 ± 0.8 keV in the south and to 13±3 keV in the north-northwest. North of the temperature peak, and slightly west of the radio contours (see Fig. 11), we find a ∼20 keV elongated sub-structure; likewise, to the south and west of the peak, we see another ∼20 keV sub-structure.

Farther out towards the outskirts, in the north-northeast direction, we see a spatially large region with a temperature of ∼10 keV. To the west of MACS J0717.5+3745, after we pass the ∼10 keV region, we observe ∼5 keV plasma surrounded by 6–7 keV plasma. Results approaching the R₂₀₀ are probably biased by foreground structures and extended emissions to the west. In the south-southeast direction, we see a filamentary structure connecting the cluster with a group of galaxies. The filament temperature is reported in Table G.1. In the map, the group itself has a temperature of 3.61 ± 0.17 keV, whereas the temperature reported in Table G.1 gives the group R₂₀₀ temperature. Based on previous lensing studies, it was determined that the filament continues farther out in this direction, even branching off into two separate filaments (Jauzac et al. 2018). Unfortunately, our X-ray data do not allow us to detect the X-ray emission associated with this part of the filamentary structure.

The measured density map traces the X-ray surface brightness image. We can also clearly see the density excess of the filament connecting the cluster and the group. Additionally, the density map appears to trace the surface brightness excess at the ‘finger’ structure to the north-northwest described by van Weeren et al. (2017). The significant temperature and pressure discontinuities that we observe in the centre of the cluster indicate the potential presence of a shock candidate.

In particular, the pressure discontinuity to the southeast of halos C and D, (see Fig. 11) indicates that we may be looking at a shock candidate. This would also be consistent with the radio emission reported by van Weeren et al. (2017). Previous discussions about the nature of these radio structures in van Weeren et al. (2017) proposed the hypothesis that these are radio relics created as a result of sub-cluster merger activity, which are seen at some peculiar projected geometry. Based on our temperature map, we estimated the Mach number for the shock candidate somewhere between ℳ = 1.7 ± 0.3 and ℳ = 2.0 ± 0.3. The study by van Weeren et al. (2017) offered an estimate of the Mach number of the shock derived from the radio spectral index of ℳ = 2.7, which is not significantly different from our X-ray values, considering the various uncertainties. Alternatively, the observed pressure peak could be due to the gravitational potential peak of the two halos (C and D). However, this scenario appears to disagree with the observed steep temperature and pressure gradients.

For a cluster in hydrostatic equilibrium, we expect to see an azimuthally symmetric and radially monotonously increasing entropy. Any departure from such a distribution is either due to the infall of subgroups hosting lower entropy material, uplift from the cluster centre due to AGN feedback or shock heating. The entropy map shows that the group in the southeast has a relatively low entropy near the centre, which gradually increases with radius. The entropy of the X-ray bright bridge and filament appears lower than that of the surrounding ambient ICM at the same radius, which would be consistent with the potential infall of stripped material. At the position of the linear radio source, we see indications of an entropy discontinuity consistent with the presence of a shock front. This linear bar region to the south of the cluster centre is also cospatial with the temperature and pressure peak and the associated X-ray shock features. Meanwhile, the characteristic ‘V’ shape, which is present in the cluster centre, seems to be associated with the colder region just north of the shock. The complex nature of the past sub-cluster merger activity and unknown geometry gives these shock discontinuities asymmetric and complex physical structures in projection.

Curious features also become apparent in the various trend-divided thermodynamical maps. In the trend-divided density map, besides the filamentary region to the south-west, three sub-clusters near the cluster centre show excess density, specifically the bar region associated with the shock front, the ‘V’ shaped emission region at the cluster centre, and the so-called ‘finger’ region to the north-west of the cluster. These regions also appear overpressured and are associated with entropy decrements in the trend-divided pressure and entropy maps. New interesting features, now visible in the trend-divided entropy map, are two high entropy regions originating from the cluster centre and going to the northeast and southwest. Looking at the radio contours present in Fig. 11, the trend of this high entropy feature near the cluster centre seems to trace the radio shock reported by van Weeren et al. (2017). However, near the outskirts of the cluster, it might be associated with the foreground structures present in the west of the cluster system.

5.2. The filament

The area of our selected rectangular filament region is 1.2 arcmin², assuming a projected geometry, we calculate a volume of an elongated rectangular cuboid inclined near the LOS based on the geometry proposed by Jauzac et al. (2012), where the assumed projected LOS depth was determined to be ∼4000 kpc.

The filament temperature of $3.1_{-0.3}^{+0.6}$ keV is lower than the cluster and the group temperatures at the same radius. The inferred filament density based on the assumed volume is (3.78 ± 0.05)×10⁻⁴ cm⁻³. The filament has a lower entropy than the cluster and group outskirts. This can also be seen in the entropy map in the bottom right of Fig. 6, where a low entropy bridge is surrounded by a gas of higher entropy.

At the cluster redshift, the critical density of the Universe is 1.7 × 10⁻²⁹ g cm⁻³. If we assume that the total baryon density is 0.044 of the critical density of the Universe and that the baryon mass fraction is 0.15 of the total mass density (Kirkman et al. 2003), we get an overdensity of our filament of 729 relative to the mean matter density of the Universe, and an overdensity of 401 relative to the critical density. Using the same volume, the baryonic mass of the filament is ∼6.1 × 10¹² M_⊙. Assuming that the baryon fraction of the filament equals the mean cosmic value, its total mass would be ∼ 7.6 × 10¹³ M_⊙.

The inferred temperature of the filament region is relatively high compared to what would be expected of a WHIM filament. The ambiguity from the projection of the filament structure provides two potential hypotheses for this. Firstly, the increased temperature might be the result of compression as the cluster and group atmospheres interact with each other. This hypothesis is supported by the fact that the calculated density for the filament region is also relatively high for what is expected from a typical WHIM filament; so if it is a true WHIM filament, it should be compressed by the merging atmospheres. Alternatively, the excess temperature could be also explained by shock heating from the WHIM accretion in the outskirts of the group atmosphere, similar to the proposed external accretion shocks (Ha et al. 2018; Gu et al. 2019). Lensing measurements of the filament structure, which reveal that it extends further to the southeast of the group, support this scenario. These types of shocks should, in theory, have high Mach numbers (ℳ ∼ 10 − 100), but have not yet been confirmed observationally since they are present in regions with very low X-ray surface brightness. Alternatively, the cooler component can be due to stripped low entropy gas from groups and galaxies seeding gas inhomogeneities in the direction of the large scale structure filament. A similar scenario was proposed for Abell 85 by Ichinohe et al. (2015).

Using the wide-field Subaru weak-lensing observations of Umetsu et al. (2014) to measure the mean surface mass density $\overline{\mathrm{\Sigma }}$ within the filamentary region, we can provide an upper limit to the total mass of the filamentary structure. We note that the mass map shown in Fig. 12 was derived using Gaussian smoothing for visualization purposes (Umetsu et al. 2014), which is not suitable for accurate mass estimation. To this end, we employ the pixelised surface mass density map Σ(θ_i) in the cluster field and its covariance matrix C_ij from Umetsu et al. (2018), derived from the combination of two-dimensional reduced shear and azimuthally averaged magnification-bias constraints using the cluster lensing mass inversion 2D (CLUMI-2D) code. Combining complementary shear and magnification measurements, this method reconstructs the underlying surface mass density field around the lens, effectively breaking the mass-sheet degeneracy (see Umetsu 2020, for a review of cluster weak lensing). The mass map Σ(θ_i) covers a field of 24 × 24 arcmin² with N_pix = 48 × 48 pixels centreed on the cluster. We find $\overline{\mathrm{\Sigma }}=(3.7\pm 1.4)\times {10}^{14}{\mathrm{M}}_{\odot }$ Mpc⁻² using an optimal (minimum variance) estimator, $\overline{\mathrm{\Sigma }}={[A^t{C}^{-1}A]}^{-1}{A}^t{C}^{-1}\mathrm{\Sigma }$ (Appendix C of Umetsu et al. 2015), where A is an N_pix × 1 mapping matrix whose elements A_i are unity for those pixels lying inside the filamentary region and zero otherwise. The uncertainty for $\overline{\mathrm{\Sigma }}$ is given by $\sigma (\overline{\mathrm{\Sigma }})={(A^t{C}^{-1}A)}^{-1/2}$. This translates into a projected mass inside the filament region of (6.8 ± 2.7) × 10¹³ M_⊙. These large uncertainties are reasonable, given that the weak lensing-based mass reconstruction has a pixel scale of 0.5 arcmin (∼190 kpc), which corresponds to only 7 pixels inside the filament region.

Fig. 12. Weak-lensing map from the CLASH survey (Umetsu et al. 2014). The white box corresponds to the same FOV as the MUSTANG2 data and the other thermodynamical maps. The green box corresponds to the zoom-in region shown in Fig. 11.

Taking these gas mass estimates (4.2 × 10¹² M_⊙) and the projected total mass measurements (6.8 ± 2.7 × 10¹³ M_⊙) at face-value, the gas mass fraction of the filament is 4–10%. This range is comparable to the upper limit of 9% determined for the filament in Abell 222/223 (Dietrich et al. 2012) and within the range of 5–10% inferred for the filament system around Abell 2744 (Eckert et al. 2015).

6. Summary

MACS J0717.5+3745 is possibly one of the most dynamically active clusters in our Universe, making it one of the most challenging galaxy clusters to study in the X-ray, owing to its geometric complexity and extremely high central temperature, combined with the many projected merging systems. Its rich radio morphology, X-ray substructures, and nearby filamentary structures have also paved the way for many interesting papers across different wavelengths.

In this paper, we present results from deep X-ray observations of the merging cluster MACS J0717.5+3745 using a novel approach to dynamically model the instrumental and astrophysical backgrounds, as well as a new method to jointly fit Chandra, XMM-Newton, and ROSAT data from the entire FOV. We applied this method to constrain the CXB and instrumental background components, before fitting the respective regions, in an attempt to push the imaging and spectroscopic capabilities of these observatories to their limit. We provided new thermodynamical maps as far out to R₂₀₀ as possible, before using joint modeling of all available data. This offered an improvement in the spatial resolution of an order of magnitude with respect to the previous spectroscopic maps. We also present a novel use of statistical model comparison methods to disentangle the complex spectral emission from several overlapping X-ray components using Bayesian analysis tools, finding a new temperature component in the filamentary structure to the S-SE of the cluster. We also present a sensitivity map of the SZ effect decrement 10″ resolution from MUSTANG-2 on the GBT. We summarise our results below.

• Thermodynamic maps were produced from Chandra, XMM-Newton, and ROSAT data of MACS J0717.5+3745, using a new method for modelling both the astrophysical and instrumental backgrounds. These maps reveal a complicated ICM structure with several subsystems. The region south-southeast of the cluster centre shows the presence of exceptionally hot 24 ± 4 keV gas and pressure discontinuities in the cluster core indicate the presence of shocks. To the south-southeast, we detected the X-ray emission of a large-scale structure filament connecting the cluster with a group.

• We described MACS J0717.5+3745 by creating a 3D cluster model assuming five spherical beta models. The geometric correction is based on the four strong lensing centres previously reported in the literature and the group at the end of the filament structure. Using this model, we calculated the effective line-of-sight length, which was then used to calculate the volume to determine the density, pressure, and entropy in the thermodynamic maps.

• Trend-divided maps reveal regions of excess density that are related to overpressured regions with complex entropy structure. These regions appear approximately co-aligned with the radio relic emission reported by van Weeren et al. (2017). The filament structure and other sub-clusters are more clearly visible after removing the cluster average trend for each physical quantity.

• The temperature peak of 24 ± 4 keV is also the pressure peak of the cluster and it is spatially offset from the SZ peaks by 0.25 arcmin, around 100 kpc, located approximately halfway between the two lensing mass peaks. We report a range for the Mach numbers for the potential shock candidate of between ℳ = 1.7 ± 0.3 and ℳ = 2.0 ± 0.3.

• We present the first MUSTANG2 SZ maps of MACS J0717.5+3745, corrected for kinetic and relativistic SZ effects using a model based on the best-fit results from Adam et al. (2017b). The SZ maps show a close agreement with the X-ray pressure peak and the kSZ model corrections are consistent with previously published works.

• Bayesian X-ray analysis methods were used to disentangle different projected spectral signatures for the filament structure, with the AIC and BIC used to select the most appropriate model to describe the various temperature components. The statistical measures indicate that the most complex model, involving CXB+cluster+group+filament components, is the best model to represent the filament structure. We report an X-ray filament temperature of $3.1_{-0.3}^{+0.6}$ keV and density 3.8 ± 0.1 × 10⁻⁴ cm⁻³, corresponding to an overdensity of the filament of 729 relative to the mean matter density of the Universe and an overdensity of 401 relative to the critical density. We estimate the baryonic mass of the filament to be ∼6.1 × 10¹² M_⊙, while its total projected weak-lensing measured mass is ∼6.8 ± 2.7 × 10¹³ M_⊙, indicating a hot baryon fraction of 4-10%.

Together, these results provide a robust thermodynamic interpretation of MACS J0717.5+3745 with a self-consistent treatment of the astrophysical and instrumental backgrounds, despite the uncertainties associated with the 3D geometry of the complex mergers. Deeper X-ray observations will push constraints further into the cluster outskirts, while joint modeling of available SZ data from all instruments together will also help break LOS projected degeneracies and refine the pressure and velocity structure.

Acknowledgments

Firstly, JPB would like to thank the reviewer for his/her comments to the different drafts of the report. The implemented changes to the paper have added much to the scientific value and potential impact of the results. JPB would also like to give a big thank you to the members of the Chandra CIAO/SHERPA help desk and the various calibration and instrument scientists for all of their help and advice. Special thanks is to be given to Nick Lee, for spending a year and over 120 emails helping me troubleshoot the development of the instrumental background fitting procedure, and for helping understand the many discovered bugs regarding the AREASCAL keywords and the background scaling in SHERPA. Special thanks is also to be given to the XMM-Newton helpdesk team and the SAS developers, for their time and patience helping address the many SAS-related problems with the development of the analysis pipeline. I would like to thank Dominique Eckert, Jelle de Plaa, and Jelle Kaastra, for several discussions and method-related motivation; Peter Boorman and Johannes Buchner for their help with using BXA; and Michal Zajaček for many helpful discussions. The scientific results reported in this article are based in part on data obtained from the Chandra Data Archive. This research has made use of software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and SHERPA; as well as NASA’s High Energy Astrophysics Software (HEASoft) packages XSPEC. The scientific results reported in this article are also based in part on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. JPB, NW, and TP acknowledge the financial support of the GAČR EXPRO grant No. 21-13491X. The material is based upon work supported by NASA under award number 80GSFC21M0002. K.U. acknowledges support from the National Science and Technology Council of Taiwan (grant NSTC 112-2112-M-001-027-MY3) and the Academia Sinica Investigator Award (grant AS-IA-112-M04). L.D.M. has been supported by the French government, through the UCA^J.E.D.I. Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.

References

Appendix A: Automatic Background Fitting

Generally, due to accumulating ionizing doses from various orbit cycles, the quantum efficiency of CCD detectors slowly degrades over time, resulting in a change in the shape of the underlying continuum due to an increase of noise in different detector channels. The shape of the Chandra background continuum has historically been fairly stable (Bartalucci et al. 2014). The shape of the XMM-Newton background continuum however, has been observed to vary between 8% and 20% between quiet periods and maximum (Katayama et al. 2004), while also being sensitive to the cleaning method associated with the cleaning and removal of background flares (Nevalainen et al. 2005). Additionally, the relative strengths of the different X-ray fluorescence line emission components also differ across the detector area due to the shielding and other nearby electronics.

For the various, relatively small, spectral extraction regions, the background model derived from the full FOV of each detector is not the most appropriate model. However, for each region on the detector, the global, underlying spectral model created from the description of the full FOV is a good first-order approximation for the start of the automatic background fitting procedure.

The background contributions are scaled proportionally to the source and background region sizes in arcmin², by the exposure times between the source and backgrounds, and also by the high energy count rates in the 10 keV to 12 keV bands. Because the instrumental backgrounds originate from the detectors and are not dependent on the effective area of the optics, the instrumental background model should not be convolved by an Ancillary Response File (ARF); however, due to peculiarities with how models are set up in Sherpa, we instead build and convolve our model with a flat-ARF, which behaves the same as no ARF when convolved. There is a standing debate with regards to which is the correct response matrix to be used for the background since, arguably, a ‘correct’ response would need to be calibrated inside the spacecraft during flight for the detector background specifically, and a separate additional response would need to be done for the soft proton (and other particle) contributions. We explored the fit results when using the photon-derived response matrix calculated from standard spectral extraction tasks and compared it with the fit results when using a diagonal-response matrix, which effectively provides infinite instrument resolution for the background, and found that the results were comparable within a few per cent.

The background fitting procedure used in this paper works iteratively, carefully thawing and fitting model parameters in sequence, and becomes more fine-tuned as each model parameter is fit for the specific detector region. The background model is first normalised with a constant before fitting the continuum emission. The continuum is then fixed before the Gaussian line emissions are then subsequently thawn and fit, adjusting their shapes and relative ratios. Finally the normalization and continuum are then refit to slightly correct for the adjusted Gaussian lines.

The automatic background fitting routine initially defaults to fitting the full instrumental background first before fitting the source spectrum. By using robust instrumental background models, the typically difficult probes into faint signatures in poor S/N and S/b data regimes (such as cluster outskirts and filamentary structures close to R₂₀₀) can be exploited by leveraging the control over the systematics in the various backgrounds, preventing a loss of information through averaging or over smoothing of subtle spectral features which might be present in the data (Ettori & Molendi 2011). This background modelling procedure is additionally motivated by a well-known bias when subtracting backgrounds using ‘wstat’ statistics that treats every bin independently from one another (See Notebook ⁴). This ‘wstat’ bias is strongly affected by how spectra are binned, especially in low-resolution data, e.g., ROSAT. Modelling the background and using ‘cstat’ decreases the overall fit uncertainties when compared to ‘wstat’ because it treats all of the bins as continuous, even in the case of low or zero-count bins.

The caveat of this background fitting method, in particular, is that the background fitting procedure begins to fail when reaching very small regions with too few background counts. In other words, when the X-ray surface brightness is high (high S/N), the regions become small, which also increases the S/b and results in poor statistics. Coincidentally, this happens to be near the size of the PSF of each instrument. In these situations, the many parameters of the complex background model cannot be properly constrained. Therefore, the model becomes inappropriate with regard to the data quality of the background data. In these cases, we simplify the complex model to adjust to the data quality. We begin to do this by first removing the instrumental Gaussian lines so that we can fit a constant and the continuum. Then we try to also remove the continuum emission, and simply fit a constant model to the background. After each of these fits are performed, the reduced statistic, Akaike Information Criterion (AIC), and Bayes Information Criterion (BIC) are computed, and the most appropriate background model is selected for the dataset.

For the very extreme cases, other background modelling methods and tools exist which could also be used; however, in this high S/b regime, a debate could be made about the need for a background model in the first place since the background can effectively be treated as noise (Buchner et al. 2014); simmonds; suzuki21.

The AIC and the BIC both give us indicators for model selection based on the computed evidence (Akaike 1974; Sugiura 1978). The small-sample size-corrected AIC and BIC are described, respectively, via

$$\begin{array}{cc}\hfill AI{C}_c& =2k-2ln(\mathcal{L})+\frac{2k^2+2k}{n-k-1},\hfill \\ \hfill BIC& =kln(n)-2ln(\mathcal{L}),\hfill \end{array}$$(A.1)

where k is the degrees of freedom in the model, n is the number of data points, and ℒ is the maximised value of the likelihood function. The non-corrected AIC equation is the same as Eq. A.1 without the final term.

As evident by the two equations, the AIC and BIC both try to account for the goodness of fit between models, with a penalty given for the number of model parameters, where respectively, the AIC penalty is 2k, versus the BIC penalty of ln(n)k. Generally speaking, this punishes model complexity to prevent favouring models which overfit the data. These values can then be used to inform model selection; however, this is only applicable when comparing nested models relative to the same fitted data. The ‘cstat’ value is approximately equal to −1/2ln(ℒ) and has been shown to be a reliable estimator for model performance even with as few as 30 counts (Cash 1979; Kaastra 2017). We can then modify Eq. A.1 to give metrics for our background model selection.

Regardless of which background is chosen, the novelty of this tool is that the forward model is built completely from first principles, is intuitive, and is easily extendable, so new missions like XRISM and NewAthena will be able to benefit from it as well. The background modelling tool is open source and will be released alongside a full pipeline for both Chandra and XMM-Newton, which handles the intermediate processing of the spectra⁵.

Appendix B: Instrumental Cross-calibration

There have been many studies exploring the cross-calibration between Chandra, XMM-Newton, and other X-ray observatories, motivated by the International Astronomical Consortium for High-Energy Calibration (IACHEC) using supernova remnants, quasars, and relaxed galaxy clusters (see Plucinsky et al. 2017, Section 6.1.2 for a detailed overview).

Using fluxes as a proxy, Nevalainen et al. (2010) found that the fluxes between Chandra ACIS and XMM-Newton pn had a large scatter; however, they were consistent within 2% in the soft band (0.5 - 2.0 keV) and that ACIS was, on average, 11% higher in the hard band (2.0 - 7.0 keV). Schellenberger et al. (2015), on the other hand, showed a difference in ACIS fluxes between 0-10% in the soft band and a 0-5% difference between ACIS and XMM-Newton MOS. A more recent study discusses an agreement on flux between ACIS and MOS below 3 keV, but the flux discrepancy rises to 10-15% above 4 keV (Madsen et al. 2017).

Regarding measured temperatures, previous studies have shown a 5-10% difference on the effective temperature between the MOS and pn detectors (Schellenberger et al. 2015). Chandra systematically has shown higher temperatures than XMM-Newton in the broad band, yet in another more recent study using Chandra and NuSTAR data, Wallbank et al. (2022) showed a 10.5% and a 15.7% temperature discrepancy between the instruments using the broad band (0.6 - 9.0 keV) and hard band (3.0 - 10.0 keV), respectively. An even more recent cross-calibration talk from the 2023 IACHEC meeting⁶ elaborated on this study further and discusses a complex relationship between the cross-calibrations between Chandra and XMM-Newton, using NuSTAR as a reference. This talk mentions that the temperature discrepancies between instruments increase with respect to measured temperatures, also noting a ∼20% flux excess for XMM-Newton in the soft band. Furthermore, there is evidence that XMM-Newton agrees better with NuSTAR for hot clusters, while Chandra agrees better for cooler clusters.

We performed several independent tests to explore the cross-calibrations between the various instruments as part of the verification of our joint fitting methodology using the instrumental background models. We found a consistent, 12-17% cross-calibration related uncertainty between Chandra and XMM-Newton between the three different observations of MACS J0717.5+3745 and a 1-3% uncertainty between the MOS and pn detectors. To do this, we fit a circular region from the high S/N regime of MACS J0717.5+3745, where there is an extremely hot (≳20 keV) plasma emission, using a collisional ionisation equilibrium (apec) plasma model multiplied by a constant. We further explored the effect that the scaling of the instrumental backgrounds has on the fit results by scaling the Chandra backgrounds by the number of counts and also by scaling by the count-rate in the 10-12 keV band. We found that the temperatures were consistent between Chandra and XMM-Newton when the background was scaled by the ratio of the observation to background high-energy particle count rates, but that the normalizations were off by around 18%. However, when the background is instead scaled by the ratio of the high-energy number of counts rather than by exposure time, the normalizations were within 5%, yet the temperatures were then off by around 22%. Despite our measured discrepancies between the different instruments, we found that our uncertainties fell within the previously reported limits from different IACHEC studies, so our new instrumental background modelling approach does not seem to be introducing any new biases to the measurements that were not already previously known or reported.

Fig. B.1 shows the statistical errors of the various instruments as a function of temperature using our joint fitting method, by extracting spectra from various regions with different temperatures. For comparison, the systematic uncertainties for temperatures of up to ∼10 keV, are estimated to be of the order of 10–15% Nevalainen et al. (2010), while for higher temperatures, based on the studies of the Bullet cluster, they are around 20% (Markevitch 2006; Wik et al. 2014).

Fig. B.1. Plot of the 1 sigma error as a function of temperature for Chandra, XMM-Newton, and the combined fitting procedure in units of keV and as percents.

Appendix C: MUSTANG2 Full FOV image

Fig. C.1. FOV of the S/N of the MUSTANG2 data shown in Fig. 11. The FOV corresponds to the white box region shown throughout the paper, while the green box, and contours are the same as described in Fig. 11.

Appendix D: MUSTANG2 Full FOV noise image

Fig. D.1. A realization of the MUSTANG2 noise map across the FOV shown in Fig. C.1. The FOV corresponds to the white box region shown throughout the paper, while the green box, and contours are the same as described in Fig. 11.

Appendix E: 3D LOS adjusted maps and projection effects - fractional differences

The magnitude and spatial distribution of the geometric correction are clearly highlighted in Fig.E.1, showing fractional differences between corrected and original maps.

Fig. E.1. Fractional differences between the corrected thermodynamic maps and the original maps, illustrating the scale of corrections due to the assumed LoS geometry.

Appendix F: 3D LOS adjusted maps and projection effects - sigma differences

Additionally, Fig.F.1 presents these differences expressed in terms of statistical significance (σ), based on the propagated uncertainties calculated from the original spectral fit for each independent bin. We caution the reader that, due to this being a highly disturbed, merging cluster with at least five sub-clusters in various stages of merger, the geometry assumed for this estimate is highly uncertain.

Fig. F.1. Differences between the corrected and original maps expressed in units of sigma (σ), calculated from the statistical uncertainties associated with each bin.

Appendix G: Best-fit spectroscopy results of the filament

All Tables

All Figures

Fig. 1. Images of the full FOV of Chandra (left) and XMM-Newton (right) of MACS J0717.5+3745. The surface brightness excess to the west/south-west are several foreground structures in between us and the cluster, indicated by the ‘foreground emission’ label. The XMM-Newton image is over-saturated in the cluster centre so that this can be seen more clearly. The green dashed circle indicates the approximate R200 of the cluster. Intensity scale is in flux per pixel (0.492 × 0.492 arcsec for Chandra, and 4.1 × 4.1 arcsec for XMM pn).

In the text

Fig. 2. Chandra (left) and XMM-Newton (right) images of MACS J0717.5+3745, restricted to the same FOV as the SZ data and thermodynamical maps, visible as the box in Fig. 1. Intensity scale is in flux per pixel (0.492 × 0.492 arcsec for Chandra, and 4.1 × 4.1 arcsec for XMM pn). Part of the XMM-Newton image is over-saturated so that the filament structure is more visible. The white lines indicate the regions used for modelling the spectra of the outskirts of the cluster and the group together with the filament. The green dashed circle indicates the approximate R200 of the cluster.

In the text

Fig. 3. Emission-free regions taken from the observations to fit the CXB components, with Chandra on the left, and XMM-Newton on the right. The box is the FOV used for Fig. 2 and the thermodynamical maps, corresponding to the total FOV of the SZ data. The green dashed circle indicates the approximate R200 of the cluster. The excluded circles are removed sources, the large excluded box region to the west corresponds to removed foreground galaxies, and another smaller box removes the emission of the filamentary structure.

In the text

Fig. 4. Maps of temperature (top left), density (top right), pressure (bottom left), and entropy (bottom right) with corresponding units of keV, cm−3×(l/1Mpc)−1/2, keVcm−3×(l/1Mpc)−1/2, and keVcm2× (l/1 Mpc)1/3. The temperature and electron density maps were created using binning with a S/N of 120 (14400 counts) and 100 (10000 counts), respectively, while the pressure and entropy maps were produced using the results from the top two S/N maps. Density maps were derived from the normalization fits assuming a 1 Mpc projected LOS depth. Spectral fits were made using data from both Chandra and XMM-Newton for the temperature, and only Chandra for the density, each using the CXB best-fit results from ROSAT, Chandra, and XMM-Newton. Fits for these maps were performed using the most appropriate background model for each region according to the AIC.

In the text

	Fig. 5. Map of the line-of-sight effective length derived from the 3D numerical simulation. The five sub-cluster centres are shown as white circles, corresponding to the different beta models shown in Table 5.
In the text

Fig. 6. Maps of temperature (top-left), density (top-right), pressure (bottom-left), and entropy (bottom-right) with corresponding units of keV, cm−3, keVcm−3, and keVcm2. The temperature and electron density maps were created using binning with a S/N of 120 (14 400 counts) and 100 (10 000 counts), respectively, while the pressure and entropy maps were produced using the results from the top two S/N maps. Density maps were derived from the normalization fits assuming an effective length determined from the best-fit results of five beta models, shown in Table 5. Spectral fits were made using data from both Chandra and XMM-Newton for the temperature, and only Chandra for the density, each using the CXB best-fit results from ROSAT, Chandra, and XMM-Newton. Fits for these maps were performed using the most appropriate background model for each region according to the AIC.

In the text

	Fig. 7. Geometry-corrected, trend-divided maps of temperature (top-left), density (top-right), pressure (bottom-left), and entropy (bottom-right).
In the text

Fig. 8. Chandra filament spectra, including all astrophysical and instrumental model components. Spectral data was averaged for each energy bin in all relevant observations, and then re-binned for better visibility. The middle and lower plots correspond, respectively, to the residuals for the astrophysical source data with respect to the full model and the instrumental background data and related instrumental background model.

In the text

Fig. 9. XMM-Newton MOS1+MOS2 filament spectra, including all model components. Spectral data was averaged for each energy bin in all relevant observations and then re-binned for better visibility. The middle and lower plots correspond, respectively, to the residuals for the astrophysical source data with respect to the full model, as well as the instrumental background data and related instrumental background model.

In the text

Fig. 10. XMM-Newton pn filament spectra, including all model components. The spectral data was averaged for each energy bin in all relevant observations, and then re-binned for better visibility. The middle and lower plots correspond, respectively, to the residuals for the astrophysical source data with respect to the full model and the instrumental background data and related instrumental background model.

In the text

Fig. 11. Comparison plots between the geometry-corrected X-ray pressure map (top left), the geometry-corrected trend-divided X-ray pressure map (top-right), followed by SZ MUSTANG2 data (bottom-left), next to a kinetic and relativistic SZ corrected version of the data, shown here in KelvinCMB units (bottom-right). The white contours are from the radio images seen in van Weeren et al. (2017) and the green contours corresponds to the MUSTANG2 data in the bottom-left. The circles represent the strong lensing sub-cluster centres given in Limousin et al. (2016), which were used to reproduce the kSZ+rSZ model used in Adam et al. (2017b), displayed in Table 7. The full FOV of the MUSTANG2 data can be seen in Appendix C.

In the text

	Fig. 12. Weak-lensing map from the CLASH survey (Umetsu et al. 2014). The white box corresponds to the same FOV as the MUSTANG2 data and the other thermodynamical maps. The green box corresponds to the zoom-in region shown in Fig. 11.
In the text

	Fig. B.1. Plot of the 1 sigma error as a function of temperature for Chandra, XMM-Newton, and the combined fitting procedure in units of keV and as percents.
In the text

	Fig. C.1. FOV of the S/N of the MUSTANG2 data shown in Fig. 11. The FOV corresponds to the white box region shown throughout the paper, while the green box, and contours are the same as described in Fig. 11.
In the text

	Fig. D.1. A realization of the MUSTANG2 noise map across the FOV shown in Fig. C.1. The FOV corresponds to the white box region shown throughout the paper, while the green box, and contours are the same as described in Fig. 11.
In the text

	Fig. E.1. Fractional differences between the corrected thermodynamic maps and the original maps, illustrating the scale of corrections due to the assumed LoS geometry.
In the text

	Fig. F.1. Differences between the corrected and original maps expressed in units of sigma (σ), calculated from the statistical uncertainties associated with each bin.
In the text