# The XMM Cluster Outskirts Project (X-COP): Thermodynamic properties of the Intracluster Medium out to $R_{200}$ in Abell 2319

V. Ghirardini, S. Ettori, D. Eckert, S. Molendi, F. Gastaldello, E.

Pointecouteau, G. Hurier, H. Bourdin

► To cite this version:

V. Ghirardini, S. Ettori, D. Eckert, S. Molendi, F. Gastaldello, et al.. The XMM Cluster Outskirts Project (X-COP): Thermodynamic properties of the Intracluster Medium out to $R_{200}$ in Abell 2319. Astronomy & Astrophysics - A&A, 2018, 614, pp.A7. 10.1051/0004-6361/201731748. hal-01815257

HAL Id: hal-01815257

https://hal.science/hal-01815257v1

Submitted on 6 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The XMM Cluster Outskirts Project (X-COP): Thermodynamic properties of the intracluster medium out to $R_{200}$ in Abell 2319

V. Ghirardini<fnref n="1" />,<fnref n="2" />, S. Ettori<fnref n="2" />,<fnref n="3" />, D. Eckert<fnref n="4" />,<fnref n="5" />, S. Molendi<fnref n="6" />, F. Gastaldello<fnref n="6" />, E. Pointecouteau<fnref n="7" />,<fnref n="8" />, G. Hurier<fnref n="9" />,<fnref n="10" />, and H. Bourdin<fnref n="11" />,<fnref n="12" />

Received 9 August 2017 / Accepted 17 January 2018

## ABSTRACT

**Aims.** We present the joint analysis of the X-ray and Sunyaev–Zel'dovich (SZ) signals in Abell 2319, the galaxy cluster with the highest signal-to-noise ratio in SZ Planck maps and that has been surveyed within our *XMM-Newton* Cluster Outskirts Project (X-COP), a very large program which aims to grasp the physical condition in 12 local ($z < 0.1$) and massive ($M_{200} > 3 \times 10^{14} M_{\odot}$) galaxy clusters out to $R_{200}$ and beyond.

**Methods.** We recover the profiles of the thermodynamic properties by the geometrical deprojection of the X-ray surface brightness, of the SZ Comptonization parameter, and accurate and robust spectroscopic measurements of the gas temperature out to 3.2 Mpc (1.6 $R_{200}$), 4 Mpc (2 $R_{200}$), and 1.6 Mpc (0.8 $R_{200}$), respectively. We resolve the clumpiness of the gas density to be below 20% over the entire observed volume. We also demonstrate that most of this clumpiness originates from the ongoing merger and can be associated with large-scale inhomogeneities (the “residual” clumpiness). We estimate the total mass through the hydrostatic equilibrium equation. This analysis is done both in azimuthally averaged radial bins and in eight independent angular sectors, enabling us to study in detail the azimuthal variance of the recovered properties.

**Results.** Given the exquisite quality of the X-ray and SZ datasets, their radial extension, and their complementarity, we constrain at $R_{200}$ the total hydrostatic mass, modelled with a Navarro–Frenk–White profile at very high precision ($M_{200} = 10.7 \pm 0.5^{\text{stat.}} \pm 0.9^{\text{syst.}} \times 10^{14} M_{\odot}$). We identify the ongoing merger and how it is affecting differently the gas properties in the resolved azimuthal sectors. We have several indications that the merger has injected a high level of non-thermal pressure in this system: the clumping free density profile is above the average profile obtained by stacking Rosat/PSPC observations; the gas mass fraction recovered using our hydrostatic mass profile exceeds the expected cosmic gas fraction beyond $R_{500}$; the pressure profile is flatter than the fit obtained by the Planck Collaboration; the entropy profile is flatter than the mean profile predicted from non-radiative simulations; the analysis in azimuthal sectors has revealed that these deviations occur in a preferred region of the cluster. All these tensions are resolved by requiring a relative support of about 40% from non-thermal to the total pressure at $R_{200}$.

**Key words.** galaxies: clusters: general – galaxies: clusters: intracluster medium – X-rays: galaxies: clusters – intergalactic medium

## 1. Introduction

Cosmic structures evolve hierarchically from the primordial density fluctuations into larger and larger systems under the action of gravity. Galaxy clusters are the largest bound structures in the Universe and the most recent products of structure formation. Baryons fall into the gravitational potential of dark matter halos and heat up to a temperature of the order of a few million Kelvin, emitting in X-rays mostly through bremsstrahlung process. In the last few years our knowledge of the physical condition of the intracluster medium (ICM) has significantly improved through the study of the Sunyaev–Zel'dovich effect (SZ; Sunyaev &

Zeldovich 1972). It arises when photons from the cosmic microwave background (CMB) are scattered by the free electrons of the ICM. The observed distortion of the CMB spectrum is directly proportional to the thermal electronic pressure integrated along the line of sight. This linear dependence implies that the SZ signal decreases more slowly than the X-ray signal, which depends quadratically on the density. The assumption that the ICM is fully thermalized and in hydrostatic equilibrium has been made in several studies (see Ettori et al. 2013, for a review). However this assumption might not be valid in cluster outskirts where the relative contribution of non-thermal to total pressure might not be negligible (e.g. Battaglia et al. 2012).

<footnote n="1">¹ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Piero Gobetti 93/2, 40129 Bologna, Italy
e-mail: vittorio.ghirardini2@unibo.it</footnote>

<footnote n="2">² INAF, Osservatorio Astronomico di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="3">³ INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="4">⁴ Max-Planck Institut für Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="5">⁵ Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland</footnote>

<footnote n="6">⁶ INAF - IASF-Milano, Via E. Bassini 15, 20133 Milano, Italy</footnote>

<footnote n="7">⁷ CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France</footnote>

<footnote n="8">⁸ Université de Toulouse, UPS-OMP, IRAP, 31400 Toulouse, France</footnote>

<footnote n="9">⁹ Centro de Estudios de Fisica del Cosmos de Aragon, Plaza San Juan 1, Planta-2, 44001 Teruel, Spain</footnote>

<footnote n="10">¹⁰ Institut d'Astrophysique Spatiale, CNRS (UMR8617) Université Paris-Sud 11, Batiment 121, 91405 Orsay, France</footnote>

<footnote n="11">¹¹ Harvard Smithsonian Centre for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="12">¹² Dipartimento di Fisica, Università degli Studi di Roma Tor Vergata, via della Ricerca Scientifica 1, 00133 Roma, Italy</footnote>
The matter distribution in the outskirts of galaxy clusters is expected to be clumpy (Nagai & Lau 2011; Vazza et al. 2013) and asymmetric (Eckert et al. 2012; Roncarelli et al. 2013), with a substantial contribution from non-thermal physics such as turbulence and bulk motion (Vazza et al. 2011), cosmic rays (Pfrommer et al. 2007), and magnetic fields (Dolag et al. 1999). Gas clumping plays an important role in the outer parts of galaxy clusters. Zhuravleva et al. (2013) showed that the density distribution inside a given shell surrounding the cluster centre can be described by a log-normal distribution modified by the presence of a high-density tail produced by the presence of clumps. It was shown that the median of this distribution coincides with the mode of the log-normal, while the mean is biased high due to the presence of clumps. Observationally Eckert et al. (2015) have confirmed this result: they conclude that the median method is able to recover the true gas density profile when inhomogeneities are present.

The XMM Cluster Outskirts Project (X-COP; Eckert et al. 2017) is a very large programme on XMM-Newton which aims to significantly increase our knowledge of the physical conditions in the outskirts of galaxy clusters. Thirteen local and massive systems have been selected on the basis of their high signal-to-noise ratio (S/N) in the Planck survey, and reported in the first catalogue (Planck Collaboration XXIX 2014).

In this paper, we focus on Abell 2319, the most significant SZ detection in the first Planck catalogue, with a S/N of 49.0 in the second Planck catalogue (Planck Collaboration XXVII 2016). Abell 2319 is a very hot and massive cluster at low redshift ($z = 0.0557$; Struble & Rood 1999). Its galaxy distribution indicates that this is a merger of two main components with a 3:1 mass ratio, the smaller system being located 10' north of the main structure (Oegerle et al. 1995). The cluster exhibits a prominent cold front SE of the main core (Ghizzardi et al. 2010) and a giant radio halo (Farnsworth et al. 2013; Storm et al. 2015).

The paper is organized as follows: in Sect. 2 we describe the reduction and analysis of X-ray data, from background modelling to spatial and spectral analysis; in Sect. 3 we present the data reduction and analysis of the Planck SZ data; in Sect. 4, we show the reconstructed profiles of the thermodynamic quantities, describe their properties, and discuss the different methods adopted to solve the hydrostatic equilibrium equation; in Sect. 5 the analysis in azimuthal sectors is illustrated; the gas mass fraction and the hydrostatic bias are shown in Sect. 6; and the summary of our main findings and our conclusions are discussed in Sect. 7.

Throughout this paper, we assume a Λ cold dark matter cosmology with $\Omega_{\Lambda} = 0.7$, $\Omega_{m} = 0.3$ and $H_{0} = 70$ km s$^{-1}$ Mpc$^{-1}$, and $E(z) = \sqrt{\Omega_{m}(1+z)^{3} + \Omega_{\Lambda}}$. At the redshift of A2319, 1 arcmin corresponds to approximately 64.9 kpc. Uncertainties are provided at the $1\sigma$ confidence level.

In the following, we refer to and plot as reference characteristic radii, $R_{500} = 1368$ kpc and $R_{200} = 2077$ kpc, that are defined at the overdensities of $\Delta = 500$ and 200, respectively, with respect to the critical value $\rho_{c} = 3H_{0}^{2} \frac{E(z)^{2}}{8\pi G}$ and using the hydrostatic mass profile (see Table 4 in Sect. 4.4).

## 2. XMM-Newton analysis

X-ray spatial and spectral analysis provides a direct probe of density and temperature of the ICM. However, the X-ray background needs to be modelled very accurately if we want to obtain accurate measurements in the outskirts of galaxy clusters where the background dominates over the signal.

### 2.1. Data reduction

The XMM-Newton Science Analysis System (XMM-SAS v15.0) and the corresponding calibration files were used to reduce the X-ray data, following the Extended Source Analysis Software analysis scheme (ESAS; Snowden et al. 2008). The presence of anomalous individual CCDs is also taken into account, removing them from the analysis. Soft proton flare periods are filtered out using the ESAS tasks *mos-filter* and *pn-filter*, thus providing clean events files. The ESAS procedure *cheese* is adopted in order to mask point sources which contaminate the field of view. Spectra, effective areas, and response files (ARF and RMF respectively) for the selected regions are extracted using the ESAS tasks *mos-spectra* and *pn-spectra*.

This procedure is applied to all seven of the observations we use in the analysis of Abell 2319: an archival central exposure, four offset observations (done specifically for the X-COP program), and two other archival exposures pointing just outside the virial radius and used to estimate the local sky background. Table 1 provides some information regarding these observations, such as the OBSID, the total and the clean exposure time, and the level of soft protons contamination obtained by comparing the measured count rate in a hard spectral band in the exposed and unexposed part of the field of view (inFOV/outFOV, Leccardi & Molendi 2008).

### 2.2. Particle background modelling

We extracted count images from the cleaned event files in the [0.7–1.2] keV energy band, where we expect to maximize the signal-to-background ratio (e.g., Ettori et al. 2010; Ettori & Molendi 2011). In Appendix A we present our method for modelling the 2D distribution and intensity of the non-X-ray background (NXB), distinguishing its different components, and computing the total NXB image in the required energy band. We briefly summarize the main steps here.

The XMM-Newton NXB is made of three separate components: the quiescent cosmic-ray induced particle background (QPB), the soft protons (SP), and a stable quiescent component (QC) whose origin is yet unknown (Salvetti et al. 2017). To model the QPB we used the unexposed corners of the EPIC cameras to estimate the QPB level in each observation. We then used filter-wheel-closed observations to model the spatial distribution of the QPB and renormalized the filter-wheel-closed data to match the count rate measured in the unexposed corners.

The residual contribution after subtraction of the QPB is split between the QC and SP components. In Appendix A we describe our method for taking the relative contribution of these two components into account. Briefly, we measure radial surface brightness profiles for a large sample of 495 blank-sky pointings and we optimize the relative contribution of these components as a function of the estimated SP contamination, imposing that the residual surface brightness profiles be consistent with a flat curve in the energy band of interest. This procedure leads to an accurate modelling of the SP contamination, as shown in Fig. A.1. The deviations from a flat profile are found to be at a level of less than 5%, thus for the remainder of the paper we adopt a systematic uncertainty of 5% of the NXB level on the measured surface brightness profile.

### 2.3. Spatial analysis

We combined the results from MOS1, MOS2, and *pn* and mosaicked all seven observations into one single image.
<table><caption>**Table 1.** A2319 observation log.</caption><thead><tr><th>Observation</th><th>OBSID</th><th>Total [ks]</th><th>MOS1 [ks]</th><th>MOS2 [ks]</th><th>pn [ks]</th><th>inFOV/outFOV</th></tr></thead><tbody><tr><td>Centre</td><td>0600040101</td><td>58.3</td><td>48.3</td><td>49.3</td><td>41.1</td><td>1.215</td></tr><tr><td>North</td><td>0744410101</td><td>36.0</td><td>23.8</td><td>24.5</td><td>19.4</td><td>1.132</td></tr><tr><td>South</td><td>0744410301</td><td>31.0</td><td>13.8</td><td>14.0</td><td>7.0</td><td>1.406</td></tr><tr><td>East</td><td>0744410401</td><td>41.9</td><td>14.4</td><td>15.4</td><td>9.5</td><td>1.346</td></tr><tr><td>West</td><td>0744410201</td><td>37.5</td><td>23.4</td><td>25.1</td><td>9.8</td><td>1.152</td></tr><tr><td>Outside</td><td>0743840201</td><td>15.0</td><td>12.1</td><td>12.3</td><td>5.7</td><td>1.261</td></tr><tr><td>Outside2</td><td>0763490301</td><td>18.0</td><td>12.9</td><td>12.8</td><td>9.0</td><td>1.253</td></tr></tbody></table>

**Notes.** The table lists: pointing name, OBSID, total exposure time, and clean exposure time for MOS1, MOS2, and *pn*, and inFOV/outFOV ratio, for the seven observations used in this work. All the observations were obtained using the medium filter, the full frame science mode for MOS, and extended full frame for *pn*.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Mosaicked and Voronoi tessellated image of A2319 in the energy band [0.7–1.2] keV corrected for the particle background. The red region is the one chosen for the estimate of the local sky background. The green circle represents the location of $R_{200}$.</figcaption></figure>

We filtered our image one more time by using the *Chandra* tool *wavdetect* in order to find the remaining point sources which contaminate the field of view but were missed by the ESAS task *cheese*. Indeed, this procedure has some difficulties finding some obvious point sources that lie near the gaps of the CCDs, or not found due to the parameters adopted.

A Voronoi tessellation algorithm (Diehl & Statler 2006) was applied on the mosaicked count image to create an adaptively binned surface-brightness map with a minimum of 20 counts per bin. The resulting Voronoi tessellated count rate map for A2319 is shown in Fig. 1.

In order to analyse spatially the cluster's image, we chose a background region located as far as possible from the cluster centre in order to have a good estimate of the sky background and to minimize the cluster contamination. We chose all the pixels in the image beyond 42 arcmin from the cluster centre as the region where we estimate the local sky background (the red region in Fig. 1). The background level is just the mean count rate in this region: $(1.82 \pm 0.06) \times 10^{-4}$ cts s$^{-1}$ arcmin$^{-2}$, in the energy band [0.7–1.2] keV (or, converting in flux using a power law spectral model with photon index $1.41: 1.46 \pm 0.05 \times 10^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$).

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Background-subtracted surface brightness profiles in the [0.7–1.2] keV energy band using the mean and median methods (red and green points, respectively). The sky background level is shown with a horizontal dashed line. The vertical dotted and dashed line represents the location of $R_{500}$ and $R_{200}$, respectively.</figcaption></figure>

The background-subtracted surface brightness profile was then computed in annular regions. We chose the annuli such that the total amount of net count rate in the [0.7–1.2] keV energy band is the same in all the regions. This choice ensures comparable statistics in all annuli. Using the ARF and RMF files for MOS2 (since the combined image was in units of MOS2), we were able to convert from count rates to fluxes. As shown in Fig. 2, we also evaluated the surface brightness from both the azimuthal mean and the azimuthal median of the brightness distribution. Following the analysis in hydrodynamical simulations on the effects of the densest substructures on the average gas density profile (Zhuravleva et al. 2013; Roncarelli et al. 2013), Eckert et al. (2016) show that the median is indeed less biased than the mean – it is a more robust estimator since it is unaffected by compact X-ray substructures filling a small fraction of the total volume – and that the ratio between mean and median can be used to estimate the relative impact of the detected clumps, providing an estimate of the level of gas clumpiness.

The electron density was then recovered using two different techniques: the “onion-peeling” technique (e.g., Ettori et al. 2010) and the multiscale technique (Eckert et al. 2016). Both assume the emission to be spherically symmetric. The latter technique also requires a super-parametric functional form for the density profile, decomposing the surface brightness in a very large number of β-models which can be individually deprojected. We obtained electron density profiles that are consistent within $0.7\sigma$, and mean relative deviation of 5% up to the virial radius (see Fig. 3).
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Top: density profile recovered from the median surface brightness profile using the multiscale and the onion peeling technique (red line and blue points, respectively). The density coming from the spectral analysis is also shown here (green points). Abundance (middle) and temperature (bottom) from the fitting of the spectra in 19 annular regions. The vertical dotted and dashed line indicates the location of $R_{500}$ and $R_{200}$, respectively.</figcaption></figure>

## 2.4. Spectral analysis

In order to recover the electron temperature and metal abundance of the X-ray emitting plasma, we performed a spectral analysis by fitting the spectra with an absorbed thermal component in the energy range [0.5–11.3] keV, and excluding the spectral regions with strong instrumental emission lines ([1.2–2.0] keV for MOS and [7.1–9.2] keV for pn), using XSPEC (Arnaud 1996).

We extracted spectra in 19 concentric annuli, defined in order to reach an approximately constant count rate in the [0.7–1.2] keV energy band. The number of net counts in the [0.5–11.3] keV energy band and the signal-to-background ratio are listed in Table B.1. We also extracted spectra from the background region indicated in Fig. 1, which is the same region used to estimate the local sky background for the spatial analysis and where there is no evidence of cluster emission.

To model spectrally the NXB component we followed Leccardi & Molendi (2008), modelling the spectra from the unexposed region of the instruments, the QPB component, using a broken power law in the energy range [0.5–11.3] keV for MOS, and in [0.5–14.0] keV for pn excluding the energy bands [7.1–9.2] keV where we observe strong instrumental emission lines.

<table><caption>**Table 2.** Components and parameters of the X-ray background adopted in our spectral analysis.</caption><thead><tr><th>Component</th><th>Normalization</th><th>Temperature</th></tr></thead><tbody><tr><td>Galactic halo</td><td>1.0 × 10<sup>-4</sup> cm<sup>-5</sup></td><td>0.35 keV</td></tr><tr><td>Local bubble</td><td>4.1 × 10<sup>-4</sup> cm<sup>-5</sup></td><td>0.15 keV</td></tr><tr><td>CXB</td><td>1.1 × 10<sup>-4</sup> photons<br>keV cm<sup>2</sup> s at 1 keV</td><td>/</td></tr></tbody></table>

We fitted the background spectra produced by the ESAS tasks *mos-back* and *pn-back*, which yields the unexposed spectrum representative of the QPB component. This fixes the parameters of the QPB background component.

Then, in the source spectrum, we restricted to a hard band above 5 keV, and we modelled the remaining particle background component using a broken power law with shape parameters (i.e., slopes and break energy) fixed according to the results obtained in other works (see Kuntz & Snowden 2008; Leccardi & Molendi 2008), leaving only normalization free. We included in the fit a thermal component – *apec* model in the X-ray spectral fitting package (XSPEC, version 12.9.1; see Arnaud 1996) with only normalization free, and using a temperature of 9.6 keV (Molendi et al. 1999), redshift of 0.0557 (Struble & Rood 1999), and 0.3 solar abundance – considering that in the hard band the emission from the cluster is small but not negligible. In this way we fixed the parameters describing the quiescent component.

We rescaled the model particle background NXB from the whole field of view to the local sky background region, leaving all the parameters of the models fixed. Normalizations were rescaled according to areas, i.e., if a spectrum comes from half of the field of view, normalizations were halved accordingly. This way the instrumental background and the contamination from soft protons are modelled in the background region. Then, we modelled the local sky background, which is the remaining source of emission in the background region. We constructed the model using three different components:

- the Cosmic X-ray Background (CXB), which is modelled as an absorbed power law with photon index fixed to 1.41 (De Luca & Molendi 2004);

- the local bubble component, which is modelled as an unabsorbed thermal model with temperature free to vary around 0.11 keV (Liu et al. 2017), redshift equal to 0, and fixed solar elemental abundance;

- the galactic halo component, which is modelled as an absorbed thermal component with temperature free to vary around 0.22 keV (McCammon et al. 2002), redshift fixed to 0, and fixed solar elemental abundance.

Using the emission model *tbabs(apec+powerlaw)+apec*, we fitted together all the spectra extracted from the background region and obtained the sky components with normalizations and temperatures listed in Table 2, which provide a flux of $1.7 \times 10^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$ in the energy band [0.7–1.2] keV.

Similarly to the background region, we extracted the spectra in the selected regions, and we rescaled the particle background model for the NXB from the field of view where they were calculated to the specific region of interest with just a change in normalization proportional to the areas. We obtained the sky background components from the background regions, rescaled according to the covered areas.

We fixed the particle and sky background. The cluster emission was modelled with a thermal component absorbed from our
<figure><img src="image_5.png" /><figcaption>Fig. 4. Comptonization map of Abell 2319 reconstructed using MILCA (Hurier et al. 2013), with an angular resolution of 7 arcmin FWHM. The black and white circles indicate the location of R$_{500}$ and R$_{200}$, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comptonization profile extracted from the SZ map. The vertical dotted and dashed line indicates the location of R$_{500}$ and 2R$_{500}$, respectively.</figcaption></figure>

own Galaxy (model tbabs · apec in XSPEC). The gas temperature, abundance, and normalization are free parameters in the spectral fit, whereas the redshift is fixed. The galactic hydrogen column density is left free to vary between 7.2 × 10$^{20}$ and 12.8 × 10$^{20}$ cm$^{-2}$, where the lower value represents the minimum of the Galactic column density due to atomic hydrogen (as tabulated in LAB HI Galactic survey in Kalberla et al. 2005) estimated over the surveyed area, and the higher value indicates the maximum column density over the same area, also corrected for molecular hydrogen as suggested in Willingale et al. (2013).

We fitted jointly all the spectra belonging to the same annulus but extracted from different observations using the C-statistics. The best-fit parameters are shown in Fig. 3 (with goodness of the fit, net counts, signal-to-background ratio, and best-fit n$_{H}$ indicated in Table B.1).

Modelling the ICM emission with a thermal component allows X-ray observations to provide a direct probe of the gas electron density, n$_{e}$. In fact, its normalization K$_{apec}$ can be written as

$$K_{\text{apec}} = \frac{10^{-14} \text{ cm}^{-5}}{4\pi D_{\text{A}}(1+z)^2} \int_V n_e n_p dV, \quad (1)$$

where the proton number density, n$_{p}$, is proportional to n$_{e}$ (n$_{p}$ ~ 0.8n$_{e}$).

We recovered the 3D profiles, temperature, and abundance by adopting the onion peeling technique (Kriss et al. 1983; Ettori et al. 2002, and references therein). Assuming a constant gas density inside each shell, we can rewrite Eq. (1) as a matrix product (using “#” to indicate it): K$_{apec}$ ∝ V #n$_{e}^2$, where V$_{i}^j$ is the geometrical volume of the jth shell intercepted by the ith annulus. By inverting this linear equation, we obtain the electron density inside each shell as n$_{e}$ ∝ √V$^{T-1}$ #K$_{apec}$. Values of the temperature and metal abundance in each shell are then obtained as Y$_{3D}$ = V$^{T-1}$ #Y$_{2D}$ · EM / V$^{T-1}$ #EM, where EM = ∫ n$_{e}^2dV is the emission measure and Y is the quantity of interest (either temperature or metallicity; for a discussion on the systematic effects

see Ameglio et al. 2007). The errors are estimated through a Monte Carlo process.

### 3. Planck analysis

The SZ effect provides a direct measurement of the thermal pressure integrated along the line of sight (Sunyaev & Zeldovich 1972). The dimensionless Comptonization parameter is defined as

$$y(r) = \frac{\sigma_{\text{T}}}{m_e c^2} \int P_e(\ell) d\ell, \quad (2)$$

where the integral is computed along the line of sight, $\ell$, at the radius, r, from the centre; $\sigma_{\text{T}}$ is the Thomson cross section, m$_{e}$ the mass of the electron, and c the speed of light.

The pressure profile is recovered from the SZ signal measured in the all-sky survey by the Planck mission (Tauber et al. 2010; Planck Collaboration I 2016). The SZ signal map is derived from the internal linear combination of the six frequency bands of the high-frequency instrument (HFI; Lamarre et al. 2010; Planck HFI Core Team 2011) on board the Planck satellite. More specifically, we made use of the Modified Internal Linear Combination Algorithm (MILCA, Hurier et al. 2013) which offers the possibility to reconstruct the targeted signal component at various scales contributed differently by the six combined input frequency maps. We therefore reconstructed a y-map for A2319 with an angular resolution of 7 arcmin FWHM (see Fig. 4).

From the y-map, we proceeded according to the method used and detailed in Planck Collaboration Int. V (2013).

We extracted the y-parameter radial profile of A2319 from our MILCA y-map, i.e., the profile was extracted on a regular grid with bins of width $\Delta\theta/\theta_{500}$ = 0.2. The local background offset was estimated from the area surrounding the cluster beyond 5 × $\theta_{500}$ = 106 arcmin. The resulting profile is shown in Fig. 5. The pressure profile was then obtained following the real space deconvolution and deprojection regularization method first described in Croston et al. (2006), assuming spherical symmetry for the cluster. The correlated errors were propagated from the covariance matrix of the y profile with a Monte Carlo procedure and led to the estimation of the covariance matrix of the pressure profile P$_{e}$ (r).

Abell 2319 is the highest signal-to-noise ratio SZ detected cluster in the Planck SZ catalogues (S/N ~ 50; see Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). Its proximity and its extension makes it fully resolved
<figure><img src="image_7.png" /><figcaption>Fig. 6. Planck correlation matrix $\rho_{X,Y}$ for the unbinned Comptonization parameter profile.</figcaption></figure>

even at the moderate angular resolution of the Planck survey, and its SZ signal extends well beyond $R_{500}$ at high significance. We thereby were able to perform an azimuthal analysis in eight azimuthally resolved sectors (see Sect. 6). The $y$ and pressure profiles in each sector were obtained as described above after masking the $y$-map and its associated error map according to the sector definition.

Due to the moderate angular resolution of the Planck survey and the oversampling implied by our sampling of the $y$-map, we introduced co-variance between the individual pixels. It cascades on the $y$ and pressure profiles computation, hence their respective covariance matrix.

In Fig. 6 we show the correlation matrix between data points, defined as

$$ \rho_{X,Y} = \frac{\Sigma(X, Y)}{\sigma_X^2 \sigma_Y^2}, $$

where $\Sigma$ indicates the covariance matrix.

Consequently, we note that points of our $y$ and SZ pressure profiles are correlated and that the respective error bars displayed in the figures of this paper only represent the square root of the diagonal of the covariance matrix. Nevertheless, when pressure is used to derive other quantities we make complete use of the whole covariance matrix, and therefore we consider any impact of the Planck PSF in our calculations.

## 4. Joint X-ray–SZ analysis of the thermodynamic properties

The profiles of the electron density estimated from X-rays and of the pressure obtained through SZ can be combined to recover all the thermodynamic quantities that define the properties of the ICM:

- the gas temperature

$$ T = P \cdot n_e^{-1}, \qquad (3) $$

- the gas entropy

$$ K = T \cdot n_e^{-2/3} = P \cdot n_e^{-5/3}; \qquad (4) $$

<figure><img src="image_8.png" /><figcaption>Fig. 7. Clumping factor radial profile for both techniques, onion peeling (blue) and multiscale (red). The pink area represents the $1\sigma$ confidence interval around the multiscale clumping factor. The black squares represent the observed value for the clumpiness in the work of Eckert et al. (2015). The vertical dotted and dashed line marks the position of $R_{500}$ and $R_{200}$, respectively.</figcaption></figure>

- the gas mass

$$ M_{\text{gas}}(<R) = 4\pi \int_{0}^{R} \rho_{g}(r') r'^2 dr', \qquad (5) $$

where the gas mass density $\rho_g = (n_e + n_p)m_u\mu$ with $m_u$ being the atomic mass unit and $\mu \approx 0.6$ the mean molecular weight in amu;

- the hydrostatic gravitating mass:

$$ M_{\text{tot}}(<r) = - \frac{r^2}{G\rho_g(r)} \frac{dP_g(r)}{dr}, \qquad (6) $$

where $G$ is the gravitational constant, and the gas pressure $P_g$ satisfies the ideal gas law $\rho_g kT/(\mu m_u) = P_g$. The gas mass fraction is then defined as $f_{\text{gas}} = M_{\text{gas}}/M_{\text{tot}}$.

### 4.1. Clumpiness profile

X-ray imaging can be directly used to estimate the level of inhomogeneities present in the ICM. The clumping factor $C = \langle n_e^2 \rangle / \langle n_e^2 \rangle$ measures the bias that affects the reconstruction of the gas density from the X-ray emission, which is directly proportional to $n_e^2$. Since we are considering the X-ray signal collected in a narrow energy range ([0.7–1.2] keV), which is almost insensitive to the gas temperature, we can directly use the results from the spatial analysis to estimate the gas clumping factor $C$.

In a first approximation, the density distribution inside a volume shell can be described by a log-normal distribution skewed by the presence of denser outliers or clumps (Zhuravleva et al. 2013; Roncarelli et al. 2013). Therefore, while the mean of this distribution tends to overestimate the gas density, the median is robust against the presence of clumps (Eckert et al. 2015), and we can estimate $C$ as the ratio of the deprojected X-ray surface brightness profiles obtained from (i) the mean of the azimuthal distribution of the counts in annuli and (ii) the median of the same distribution. The resulting profile is shown in Fig. 7 and indicates a $\sqrt{C}$ of about 1.1 at $R_{200}$.

However, we can only detect clumps that are resolved by XMM-Newton, i.e. clumps on scales larger than the PSF half energy width (~17 arcsec $\approx 18.4$ kpc, for MOS$^1$; see also Read

<footnote>$^1$ https://heasarc.nasa.gov/docs/xmm/uhb/onaxisxraypsf.html</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Two-dimensional temperature profiles using X-ray spectral data (blue points; thick error bars represent the systematic uncertainty as estimated in see Sect. “Systematic uncertainties on the temperature profile,” and thin error bars indicate the total uncertainties); the pressure from SZ divided by density from X-ray projected on the plane of the sky (red points); and the projection of the reconstructed temperature from the backward technique, which makes use of both X-ray and SZ data, on an NFW mass model (black line). The grey shaded area is the $1\sigma$ confidence region around the backward result. The vertical dotted and dashed line marks the position of $R_{500}$ and $R_{200}$, respectively.</figcaption></figure>

et al. (2011)). This implies that clumped structures below this scale might still bias our measured thermodynamic quantities.

## 4.2. Temperature profile

Similarly to what was done for the pressure, we were able to recover the ICM temperature profiles in two different ways: (i) from the spectral analysis ($T_X$), as detailed in Sect. 2.4 and (ii) by dividing $P_{SZ}$ by the gas density $n_e$ recovered from the deprojection of the X-ray surface brightness ($T_{SZ}$). These values can be compared with the profile $T_{NFW}$ that is obtained from the best-fit mass model (see Sect. 4.4) by requiring that the hydrostatic equilibrium holds between the cluster potential and the observed gas density profile. We note that $T_{NFW}$ is not independent of the other two profiles because the best-fit mass model is obtained by fitting both the measured $T_X$ and $P_{SZ}$. In order to obtain a meaningful comparison with $T_X$, we compute a spectroscopic projection (see Mazzotta et al. 2004; Morandi et al. 2007) of the 3D quantities $T_{SZ}$ and $T_{NFW}$. The good agreement among these profiles is shown in Fig. 8.

Because the pressure gradient in the first point is washed out from the Planck’s beam of about 7 arcmin, we note that the pressure in this point is underestimated, and therefore the temperature $T_{SZ}$ is also underestimated with respect to $T_X$.

Systematic uncertainties on the temperature profile. We constrain the projected spectroscopic temperature (see Sect. 2.4) with a relative statistical uncertainty ranging between 1% and 6% (median value: 2%). It is thus critical to evaluate the role of possible systematics in our measurements. In order to calculate some of the most relevant systematic uncertainties affecting our temperature measurements, we re-estimate the spectral temperature using several different methods. Our reference temperature measurement is the one calculated using both MOS and pn data, leaving $n_H$ free to vary within a defined narrow range, and fixing the parameters of the background model. By changing all these quantities one by one, we estimate the level of systematic errors that affect our measurements. In detail, we calculate the spectral temperatures in four ways: (i) only using counts collected from MOS, (ii) only using counts from pn, (iii) fixing $n_H$ to the LAB value (Kalberla et al. 2005), and (iv) allowing

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Temperature and abundance profile adopting different techniques; this provides an estimate of the systematic error affecting our measurement.</figcaption></figure>

the background parameters (normalizations) to vary within ±5% of the best-fit values. We show in Fig. 9 the results of this procedure. Finally, at each radial point, we estimate the systematic error using the standard deviation of the values measured with all the different methods. This error is then added, in quadrature, to the statistical error and propagated through the entire analysis. The relative systematic error ranges between 1.4% and 9.1%, apart from the outermost radial point where we measure a value of 19%.

## 4.3. Pressure profile

If the galaxy cluster is not affected by an ongoing merger generating shocks through the ICM, the pressure is the thermodynamic quantity that presents a smoother spatial distribution along the azimuth. It is described by a “universal” form (Nagai et al. 2007; Arnaud et al. 2010),

$$ \frac{P(x)}{P_{500}} = \frac{P_0}{(c_{500}x)^\gamma [1 + (c_{500}x)^\alpha]^{\frac{\beta-\gamma}{\alpha}}} \quad (7) $$

where

$$ P_{500} = 1.65 \times 10^{-3} \text{ keV cm}^{-3} \left( \frac{M_{500}}{3 \times 10^{14} M_{\odot}} \right)^{2/3} E(z)^{8/3} \quad (8) $$

and $x = R/R_{500}$; $\gamma$, $\alpha$, and $\beta$ are the central slope, the intermediate slope, and the outer slope defined by a scale parameter $r_s = R_{500}/c_{500}$ ($R \ll r_s$, $R \sim r_s$ and $R \gg r_s$ respectively); and $P_0$ is the normalization. The values of $R_{500}$ and $M_{500}$ adopted here are presented in Table 4 (see Sect. 4.4). We list in Table 3 our best-fit values, using the entire available radial range to find the best fit.
<table><caption>**Table 3.** Best-fit parameters of the pressure profile using the functional form introduced by Nagai et al. (2007).</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>γ</th><th>α</th><th>β</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>Arnaud et al. (2010)</td><td>8.40</td><td>1.18</td><td>0.31</td><td>1.05</td><td>5.49</td><td>—</td><td>—</td></tr><tr><td>Planck Collaboration Int. V (2013)</td><td>6.41</td><td>1.81</td><td>0.31</td><td>1.33</td><td>4.13</td><td>—</td><td>—</td></tr><tr><td>SZ+X</td><td>7.7 ± 2.0</td><td>1.34 ± 0.22</td><td>0.47 ± 0.07</td><td>1.05</td><td>3.80 ± 0.22</td><td>2.62</td><td>69</td></tr><tr><td>SZ</td><td>9.6 ± 5.8</td><td>1.10 ± 0.35</td><td>0.23 ± 0.23</td><td>1.05</td><td>4.50 ± 0.47</td><td>3.47</td><td>9</td></tr></tbody></table>

**Notes.** “SZ+X” refers to the best fit done on the best-fit mass model pressure profile (see Sect. 4.4), while “SZ” refers to the best fit done only on the $P_{SZ}$.

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Rescaled pressure profile in units of $R_{500}$. The grey dashed line represents the universal pressure profile. The blue and the pink lines represents the best fits using the functional form introduced by Nagai et al. (2007) done on $P_{NFW}$ and $P_{SZ}$, respectively. The dotted and dashed vertical line represents the position of $R_{500}$ and $R_{200}$, respectively. In the bottom panel we show the ratio of $P_{SZ}$, $P_X$, and $P_{NFW}$ with the universal pressure profile (Arnaud et al. 2010).</figcaption></figure>

The electronic pressure can be directly recovered from the Comptonization profile (see Eq. (2); $P_{SZ}$), and from the deprojection of X-ray measurements of the temperature and density profiles of the emitting electrons ($P_X$). We can also estimate the pressure profile required from the best-fit mass model to satisfy the hydrostatic equilibrium ($P_{NFW}$, see Sect. 4.4). As we show in Fig. 10, these 3D pressure profiles agree well within their statistical errors.

We rescale the pressure profile by $P_{500}$ and fit it with the universal functional form (Nagai et al. 2007). The best-fitting results are listed in Table 3. The comparison with the results of Planck (Planck Collaboration Int. V 2013) and Arnaud et al. (2010) is shown in Fig. 10. We observe that the pressure profile in A2319 is well above the other two profiles, in particular in the outskirts, with values higher by about a factor of ~3.5 at $R_{200}$, which is ~2σ away from the Planck envelope (Planck Collaboration Int. V 2013).

We have also adopted a new technique (Bourdin et al. 2017) in order to evaluate the impact of the anisotropies in the Compton parameter detected in the outskirts of A2319 on the reconstructed pressure profile, and conclude that these anisotropies cannot explain the observed excess.

### 4.4. Hydrostatic mass

The total mass profile of the cluster is reconstructed by solving the hydrostatic equilibrium Eq. (6) (HEE, Binney & Tremaine 2008). In this work, we use three different methods to solve this equation and recover the hydrostatic mass profile

(e.g., Ettori et al. 2013): the backward method, the forward method, and a non-parametric method.

The backward method follows the approach described in Ettori et al. (2010, 2017) and, assuming a mass model with few free parameters (generally two), minimizes a likelihood function by comparing the predicted and observed profiles of some interesting physical quantities (e.g., temperature) to constrain these parameters. In the present analysis, we assume a Navarro–Frenk–White profile (NFW, Navarro et al. 1997) for the total mass (a more extensive discussion on the best-fitting mass models will be presented in a forthcoming publication), and constrain its two parameters, concentration and scale radius (or $R_{200}$), using both the projected temperature profile from X-ray spectral analysis and the thermal pressure profile from the SZ analysis, and maximizing the likelihood described in Appendix D.

In Fig. 11, we show the best-fit results obtained using this method to constrain the parameters of the mass model, using the median method and the multiscale technique to obtain the density profile. Very consistent results are obtained by adopting different methods to recover the input profiles of the gas temperature and density (see Table 4). We indicate with the subscript “NFW” the thermodynamic quantities corresponding to the best-fit mass model.

In the forward method, functional forms are used to fit the thermodynamic quantities, density, pressure, and temperature. Then, HEE is directly applied in order to compute the total mass radial distribution. Errors are estimated through a Monte Carlo process. As mentioned in Sect. 2.3, we use the multiscale approach (Eckert et al. 2015) to fit the emissivity profile which yields directly the fitted density functional form. We use a six-parameter functional form (Vikhlinin et al. 2006) to fit the temperature, and a five-parameter generalized NFW (Nagai et al. 2007) for the pressure. We combine in several ways the profiles of the thermodynamic quantities (density, pressure, and temperature), as detailed in Table 4, making use of a joint likelihood (see Appendix D) when all three quantities are fitted together. It is worth noticing that, while measurements of the gas density and pressure are available up to ~$R_{200}$, direct spectral estimates of the temperature are limited to regions below $R_{500}$, defining the radial range where the mass profile is more reliable in this case.

Due to the good quality data both from X-rays and SZ, we can also implement a non-parametric method in order to recover the total mass profile. We just insert pressure and density in the HEE, and we calculate the pressure derivative using a three-point quadratic Lagrangian interpolation. We note that the errors relative to this method are represented by a covariance matrix since we are using the SZ pressure profile, and therefore what is shown as an error bar in the plot is just the square root of the diagonal terms.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Left: contour plot with confidence regions at 1, 2, and 3σ (solid lines) applying the backward approach to solve HEE in order to constrain the parameters of the NFW mass model; using as inputs the multiscale technique on the median emissivity profile to obtain the density, the pressure from the direct deprojection of the y-parameter radial profile, and the temperature from the spectral analysis. Right: gas mass and total mass profile recovered using the backward approach (blue and red curves, respectively). The black crosses represent the total mass profile obtained using a non-parametric method and the green cross that obtained by applying the forward method on temperature and density profiles. The dotted and dashed vertical line marks the position of R<sub>500</sub> and R<sub>200</sub>, respectively.</figcaption></figure>

<table><caption>**Table 4.** Best-fitting results on the mass model using the techniques specified in the first column.</caption><thead><tr><th>Technique</th><th>Data</th><th><i>M</i><sub>200</sub> (10<sup>14</sup> <i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>200</sub> (kpc)</th><th><i>M</i><sub>500</sub> (10<sup>14</sup> <i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>500</sub> (kpc)</th></tr></thead><tbody><tr><td><strong>Backward</strong></td><td><strong><i>P</i><sub>SZ,NO3</sub> + <i>T</i><sub>X</sub> – Median</strong></td><td><strong>10.7 ± 0.5</strong></td><td><strong>2077 ± 33</strong></td><td><strong>7.7 ± 0.4</strong></td><td><strong>1368 ± 17</strong></td></tr><tr><td>Backward</td><td><i>P</i><sub>SZ</sub> + <i>T</i><sub>X</sub> – Median</td><td>10.6 ± 0.5</td><td>2071 ± 32</td><td>7.5 ± 0.3</td><td>1357 ± 13</td></tr><tr><td>Backward</td><td><i>P</i><sub>SZ,NO3</sub> + <i>T</i><sub>X,SYS</sub> – Median</td><td>10.3 ± 0.7</td><td>2047 ± 47</td><td>7.4 ± 0.4</td><td>1350 ± 24</td></tr><tr><td>Backward</td><td><i>P</i><sub>SZ</sub> + <i>T</i><sub>X,SYS</sub> – Median</td><td>10.5 ± 0.5</td><td>2062 ± 34</td><td>7.3 ± 0.3</td><td>1347 ± 18</td></tr><tr><td>Forward</td><td><i>P</i><sub>SZ</sub> only – Median</td><td>9.4 ± 0.5</td><td>1984 ± 40</td><td>7.4 ± 0.4</td><td>1353 ± 25</td></tr><tr><td>Forward</td><td><i>T</i><sub>X</sub> only – Median</td><td>/</td><td>/</td><td>7.3 ± 0.1</td><td>1343 ± 5</td></tr><tr><td>Forward</td><td><i>P</i><sub>SZ</sub> + <i>T</i><sub>X</sub> – Median</td><td>8.3 ± 0.3</td><td>1906 ± 20</td><td>7.8 ± 0.2</td><td>1375 ± 11</td></tr><tr><td>Forward</td><td><i>P</i><sub>SZ</sub> + <i>T</i><sub>X</sub>, β fixed – Median</td><td>8.5 ± 0.6</td><td>1923 ± 48</td><td>7.7 ± 0.4</td><td>1368 ± 26</td></tr><tr><td>Forward</td><td><i>P</i><sub>SZ,NO3</sub> + <i>T</i><sub>X,SYS</sub>, β fixed – Median</td><td>7.7 ± 0.7</td><td>1859 ± 59</td><td>7.4 ± 0.6</td><td>1354 ± 37</td></tr><tr><td>Forward</td><td><i>P</i><sub>SZ,NO3,SYS</sub> + <i>T</i><sub>X</sub> – Median</td><td>8.3 ± 0.3</td><td>1907 ± 26</td><td>7.8 ± 0.3</td><td>1373 ± 18</td></tr><tr><td>Non-parametric</td><td><i>P</i><sub>SZ</sub> – Median</td><td>9.3 ± 1.1</td><td>1979 ± 78</td><td>6.7 ± 0.5</td><td>1307 ± 33</td></tr><tr><td>Backward</td><td><i>P</i><sub>SZ</sub> + <i>T</i><sub>X</sub> – Mean</td><td>10.2 ± 0.5</td><td>2040 ± 35</td><td>7.3 ± 0.3</td><td>1346 ± 17</td></tr></tbody></table>

**Notes.** In the second column the data used to constrain the mass are listed; $P_{SZ}$ and $T_X$ refer to the SZ pressure and the X-ray temperature, respectively; the subscript “NO3” indicates that the first three Planck points were not used in the analysis; the subscript “SYS” indicates that the systematic uncertainties on the X-ray temperature are added in quadrature to the statistical errors in evaluating the $\chi^2$ (see Sect. “Systematic uncertainties on the temperature profile”); “Median” or “Mean” refers to how we computed the X-ray emissivity; “β fixed” indicates that the outer slope of the pressure profile is fixed to the best-fit value of the Planck Collaboration. In the other four columns, we quote the results on $M_{200}$, $R_{200}$, $M_{500}$, and $R_{500}$, respectively. In the first row, we indicate our reference values in bold. The last two rows present the mass reconstructed using the mean density profile, and propagating the statistical error on the temperature profile only (see Sect. “Systematic uncertainties on the temperature profile”). $R_\Delta$ is defined as $(\frac{M(R)}{4/3\pi\rho_{\Delta}})^{1/3}$.

The recovered mass profiles are shown in Fig. 11. They are all compatible within their respective error bars at the characteristic overdensities of 500 and 200.

*Systematic uncertainties on the hydrostatic mass.* In Table 4, only the statistical error on $M_{200}$ is quoted (with a relative uncertainty of about 4.7%). In this section, we evaluate the impact of some of the systematic uncertainties that affect the mass reconstruction.

The ability of the particle background model to reproduce a flat surface brightness profile when applied on blank field observations is a source of systematic uncertainty caused by the adopted procedure. As we discussed in Sect. 2, adopting

the background model described in Appendix A, we are able to reduce the systematic deviation from a flat profile below 5%. We account for this by adding 5% of the background level as an extra error in the surface brightness profile.

The results obtained by applying different methods and techniques are shown in Table 4. We estimate the level of the systematic uncertainties on the mass measurement at $R_{500}$ and $R_{200}$ of about 3.9% and 8.4%, respectively, by measuring the relative scatter around the reference value.

Another source of systematic uncertainty comes from the choice of the background region, defined in an area concentrated to the west of the cluster. Considering that A2319 has
an angular extension of ~1°, cosmic variance can influence the analysis, especially in the outskirts. Using the absorbed thermal model tbabs(apec), and fixing the parameters of the apec component, we vary the hydrogen column density alone by adopting the values of $n_H$ in regions located to the north, west, east, and south, as far as possible from the centre (at a distance of 33, 55, 36, and 39 arcmin, respectively) and remeasure the conversion factor between the count rate and the surface brightness maps. This procedure allows us to measure a relative deviation of 2% on the surface brightness, which translates into an effect of about 1.4% on the gas density and 1.1% on the mass measurement.

We therefore estimate that the total systematic uncertainties are at the level of 4.18% and 8.5% at $R_{500}$ and $R_{200}$, respectively, implying that the reference values for the hydrostatic mass are, at $R_{500}$ and $R_{200}$, respectively:

$$M_{500} = 7.7 \pm 0.4^{\text{stat.}} \pm 0.3^{\text{syst.}} \times 10^{14} M_{\odot},$$

$$M_{200} = 10.7 \pm 0.5^{\text{stat.}} \pm 0.9^{\text{syst.}} \times 10^{14} M_{\odot}.$$

### 4.5. Entropy profile

The entropy profile is recovered through the gas pressure and temperature profiles via Eq. (4). Entropy is a fundamental quantity that can track the thermal history of a cluster: it always rises when a heat flow occurs, and in the presence of only non-radiative processes it is expected to follow a power law with characteristic slope of 1.1 (Tozzi & Norman 2001; Voit et al. 2005). Deviations from this power law are observed in the central regions, requiring an entropy “floor” within ~100 kpc that is expressed through the formula (Cavagnolo et al. 2009)

$$K = k_0 + k_{100} \left( \frac{R}{100 \text{ kpc}} \right)^{\alpha} \quad (9)$$

The central entropy ($k_0$) measured with the fit in Eq. (9) is 75 ± 13 keV cm² (see Table 5), suggesting that A2319 does not possess a relaxed, cool core (e.g. Cavagnolo et al. 2009, define a CC when $k_0 < 50$ keV cm²).

However, non-radiative simulations show that the self-similar behaviour is reproduced only once entropy is rescaled by a proper quantity defined with respect to the critical density (Voit et al. 2005)

$$K_{500} = 106 \text{ keV cm}^2 \left( \frac{M_{500}}{10^{14} M_{\odot}} \right)^{2/3} E(z)^{-2/3} f_{\text{b}}^{-2/3}, \quad (10)$$

where $f_{\text{b}} = 0.15$ is the universal baryon fraction. Non-radiative simulations (Voit et al. 2005) predicts that the power law describing the entropy profile is:

$$\frac{K(R)}{K_{500}} = 1.42 \left( \frac{R}{R_{500}} \right)^{1.1} \quad (11)$$

In order to accommodate the flattening of the entropy profile observed in many disturbed systems, we add a constant to a simple power law:

$$\frac{K(R)}{K_{500}} = k_0 + k_{500} \left( \frac{R}{R_{500}} \right)^{\alpha} \quad (12)$$

In Fig. 12, we plot the measured entropy profiles, also rescaled according to Eq. (12). In Table 5, we show the best-fit results on the data using Eqs. (9) and (12). We observe that the entropy profile has a shallower slope with respect to the value predicted by the simulations (Voit et al. 2005).

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Top: entropy profiles obtained from the three different methods described in Sect. 4.5. The dashed magenta line represents the best fit obtained on the $K_{\text{NFW}}$ data using Eq. (9). Bottom: entropy profiles rescaled by $K_{500}$. The dashed magenta line represents the best fit obtained on the $K_{\text{NFW}}$ data using Eq. (12). The green lines represent the prediction from Voit et al. (2005). The dashed pink lines are the best fit using Eqs. (9) and (12). The vertical dotted and dashed line represents the location of $R_{500}$ and $R_{200}$, respectively.</figcaption></figure>

<table><caption>**Table 5.** Best-fit results for the model of the entropy profile using the three different rescalings described in Sect. 4.5.</caption><thead><tr><th></th><th>k<sub>0</sub></th><th>k<sub>100/500</sub></th><th>α</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>Eq. (9)</td><td>75 ± 13</td><td>190 ± 12</td><td>0.82 ± 0.03</td><td>129</td><td>70</td></tr><tr><td>Eq. (12)</td><td>0.055 ± 0.010</td><td>1.17 ± 0.02</td><td>0.82 ± 0.03</td><td>124</td><td>70</td></tr></tbody></table>

## 5. Analysis in azimuthal sectors

### 5.1. Thermodynamic properties

Considering the high signal-to-noise ratio of our X-ray and SZ datasets, we can perform the analysis presented in the previous sections in each of the eight azimuthal sectors with width of 45° that we define in Fig. 13. The analysis performed in the sectors allows us to measure the azimuthal variance of the physical quantities and to assess which cluster regions are more relaxed. Indeed, by dividing the observed count rate map in Fig. 1 with a cluster model with perfect spherical symmetry and emission equal to the azimuthal median surface brightness profile, we can identify where an excess in the emission due to the ongoing merger is located. As shown in Fig. 13, this excess is concentrated in the NW region (sectors 1, 2, and 3, in particular).

We show the profiles of the thermodynamic properties recovered in eight angular sectors in Fig. 14.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Left: same as Fig. 1. The white sectors represents the eight regions analysed separately, each marked by a identification number. Right: residual image obtained by dividing the flux image by the model image reconstructed from the median method. The small white circle represents the centre of the cluster, and the big white circle represents the position of R$_{200}$. The red sector represents the region which shows a clear excess in the residual map.</figcaption></figure>

In the X-ray surface brightness, we identify various features specific to each sector:

- Sector 1 has an excess in emission starting above 200 kpc with a small radial extent of about 100 kpc. This excess is due to a contamination of the merging component in this cluster, located 10 arcmin NW.

- Sector 2 also has a significant excess in the X-ray emission. This excess is located in the region where Oegerle et al. (1995) found the merging component in A2319, and has a radial extent that is quite large, from 200 to 800 kpc.

- Sector 3 has an emission slightly higher than the azimuthal average up to 1 Mpc, where a sharp transition is present reconciling the surface brightness with the azimuthally averaged value. This sector shows evidence of a non-negligible contamination from the merger.

- Sectors 4 and 5 are quite regular, with a behaviour very similar to the azimuthally averaged profile.

- Sector 6 shows the cold front that was detected in Ghizzardi et al. (2010) and located in the SE region, about 200 kpc ≈ 3 arcmin from the cluster centre.

- Sectors 7 and 8 are the most regular, and reproduce very well the combined surface brightness profile.

The pressure profile obtained from the deprojected SZ signal in each sector (see Fig. 14) shows clearly that this is the quantity least affected by the dynamical history of the cluster. For instance, the merging event (Oegerle et al. 1995) happening in the NW (sector 2) with mass ratio 3:1 is well resolved in the surface brightness/density profile, but it is not evident in the pressure profile (sector 3 has the highest values in the pressure profile, while sectors 1 and 2 are slightly below the azimuthally average profile), suggesting that the merger induced some shocks that have already propagated through the ICM and are at least partially thermalized, inducing a reasonably small scatter in the pressure profile at R$_{200}$ (see Fig. 14).

From the spectral analysis, we observe in sector 2 that the gas temperature reaches values below those measured in the

azimuthally averaged profile between 300 and 800 kpc. In sectors 1 and 3, the temperature behaves similarly, but over a narrower radial range. These radial variations can be explained by a low-temperature component contaminating sectors 1, 2, and 3 at intermediate radii. This can be associated with the accreting substructure visible in the residual map (Fig. 13), which is merging with the main cluster halo. Over the same region, corresponding to the merging component at about 500 kpc in sector 2, we also observe an increase in the metal abundance correlated to the gas at the lower temperature.

In Fig. 14, we show the entropy profiles obtained by solving the HEE with the backward method (a comparison between the profiles estimated with different methods is shown in Fig. E.3). The entropy measured in sector 2 is well below the mean value estimated in the cluster, while sector 1 and 3 are just slightly below the mean value. This suggests that a substructure with a low-entropy gas is still accreting into the cluster's halo as residual of the ongoing merger.

## 5.2. Azimuthal scatter and clumpiness

The azimuthal scatter of the recovered thermodynamic quantities is defined at each radius r as

$$ \sigma_Q(r) = \sqrt{\frac{1}{N} \sum_{i=1}^{N} \left( \frac{Q_i(r) - \bar{Q}_i(r)}{\bar{Q}_i(r)} \right)^2}, \quad (13) $$

with $Q = \{n, P, T, K, M_{tot}, M_g, f_g\}$. The profiles of the azimuthal scatter are shown in Fig. 15.

As a general trend, we expect that $\sigma_Q(r)$ should increase monotonically with radius because as they move outward the considered radial points should be less virialized. Although this is generally observed, some other features also appear. For instance, at intermediate radii (~600 kpc) there is a clear increase coincidentally with the cluster centric location where the merger is taking place. Moreover, there is a particular radial location between R$_{500}$ and R$_{200}$ where the azimuthal scatter reaches a minimum. This point suggests the radial
<figure><img src="image_15.png" /><figcaption>Fig. 14. Top left: surface brightness profiles for the eight sectors using the median method. The thick black line is obtained from the whole image analysis. The dotted and dashed vertical line marks the position R$_{500}$ and R$_{200}$, respectively. Top right: SZ pressure profiles for the eight sectors overplotted on top of the pressure profile for the whole cluster (black line). (Centre) Two-dimensional temperature (left) and abundance (right) profiles for the eight sectors analysed. Bottom left: reconstructed entropy profiles for all sectors using the backward method. Bottom right: gas fraction profiles recovered applying the backward technique. The thick black line is the result for the azimuthally averaged profile. The dotted and dashed vertical line marks the position R$_{500}$ and R$_{200}$, respectively.</figcaption></figure>

extension of the influence of the merger on the thermodynamic quantities.

Using this information, we can improve the characterization of the properties of the observed clumpiness in the gas density. As described in Roncarelli et al. (2013), the clumping factor of the gas (see Sect. 4.1) is expected to have two major contributors: some individual clumps and large-scale accretion patterns. The latter is described by the residual clumping C$_{R}$, which, following Roncarelli et al. (2013), can be estimated as

$$C_R(r) = 1 + \frac{\sigma}{\sigma_0} + \frac{r}{r_0}, \qquad (14)$$

where r = R/R$_{200}$; σ is the azimuthal scatter of the density n or of the Comptonization parameter y; and σ$_{0}$ and r$_{0}$ are estimated from simulations (Roncarelli et al. 2013)

$$-(\sigma_0, r_0) = (16.02, 5.87) \text{ for X-ray density,}$$

$$-(\sigma_0, r_0) = (2.83, 8.25) \text{ for SZ Comptonization parameter.}$$

We compare the estimated clumpiness with the residual clumpiness C$_{R}$ in Fig. 15. We observe that the measured clumping factor, both X-ray and SZ, only slightly exceeds the estimated C$_{R}$ over the entire radial range, suggesting that large-scale asymmetries account for most of the clumpiness measured.

Moreover, the clumpiness profile in Fig. 7 shows a clear excess at intermediate radii. We interpret this excess as the presence of the merger component in the NW direction. We evaluate again the clumpiness, after masking out sectors 1, 2, and 3, which are more affected by the presence of the merger. As we show in Fig. 15, the excess in the clumping factor at intermediate radii disappears and the total clumpiness at R$_{200}$ decreases to 1.05.
<figure><img src="image_16.png" /><figcaption>Fig. 15. Top: azimuthal scatter in the thermodynamic profiles: gas density and gas mass profiles are obtained from the X-ray spatial analysis; the pressure profile is the result of SZ data analysis; gas entropy and temperature are obtained by combining SZ pressure and X-ray density; the total mass is reconstructed by solving the hydrostatic equilibrium equation using the forward approach. The vertical dotted and dashed line represents the location of R$_{500}$ and R$_{200}$, respectively. Centre left: total measured clumpiness (see Sect. 5.2; black line, shaded region represents 1σ uncertainty) compared with the estimated residual clumpiness using X-ray density (blue line) and the SZ Comptonization parameter (green line). Centre right: same as Fig. 7, but removing the problematic sectors (1, 2, and 3) from the analysis. The features present in the whole clumpiness profile disappear almost completely. The dotted and dashed vertical line marks the position of R$_{500}$ and R$_{200}$, respectively. Bottom left: total measured clumpiness (see Sect. 5.2; black line, shaded region represents 1σ uncertainty) compared with the estimated residual clumpiness using X-ray density (blue line) and SZ Comptonization parameter (green line). Bottom right: same as Fig. 7, but after removing the merging region in the problematic sectors (1, 2, and 3) from the analysis. The features present in the whole clumpiness profile disappear almost completely. The dotted and dashed vertical line marks the position of R$_{500}$ and R$_{200}$, respectively.</figcaption></figure>

# 6. Characterizing the hydrostatic bias

## 6.1. Gas mass fraction and the non-thermal contribution

Since galaxy clusters originate from large regions of the primordial Universe, their baryon fraction is expected to be close to the universal fraction.

The gas mass fraction, $f_g = M_g/M_{tot}$, in massive galaxy clusters represents most of the baryons accreted in the dark matter halo and is a good proxy of the cosmic baryonic budget, which enables us to use galaxy clusters as a cosmological probe (e.g., Ettori et al. 2002, 2009),

$$ \frac{\Omega_b}{\Omega_m} \cdot b = f_g + f_{star}, \qquad (15) $$

where $\Omega_b$ and $\Omega_m$ are the cosmological baryon and matter density, b is the depletion factor that accounts for the cosmic baryons which thermalize in the cluster's potential, and $f_{star}$ is the stellar mass fraction. Here, we adopt the cosmological parameters $\Omega_b = 0.045$ and $\Omega_m = 0.3089$ estimated from Planck Collaboration XIII (2016); we assume from numerical simulations $b = 0.85$ and $0.87$ (with a standard deviation of 0.03) at R$_{500}$ and R$_{200}$, respectively (e.g., Planelles et al. 2013); and consider $M_{star}/M_{gas} = 0.069$ from optical measurements in

nearby systems (Gonzalez et al. 2013). We predict, thus, a gas mass fraction $f_g$ of 0.125 and 0.128 at R$_{500}$ and R$_{200}$, respectively.

However, we measure a gas fraction, already corrected for the resolved gas clumpiness using the median profile, that reaches values well above the expected $f_g$ at $r > R_{500}$ (see Figs. 14 and 15). We argue for the role of the non-thermal pressure contribution to the estimate of the total mass in lowering the measured gas fraction.

Indeed, Abell 2319 is in a merging state (Oegerle et al. 1995), with the presence of a giant radio halo (Farnsworth et al. 2013; Storm et al. 2015) that supports this scenario. The measured gas fraction can then be biased high as a consequence of the phenomena (like gas turbulence and bulk motion) that occur during a merger and that are not accounted for in the calculation of the hydrostatic mass, causing an underestimate of the halo mass.

Before proceeding in quantifying the amount of non-thermal pressure support, we note (from the analysis in azimuthal sectors) that the substructure that is merging with the main halo is also able to disturb the system on a much larger scale by enhancing the measured surface brightness up to ~1 Mpc. The net effect is to increase the gas mass by about 10% and so the relative amount of non-thermal pressure in the outskirts. To obtain an estimate of the contribution of the non-thermal pressure unbiased from any evident merger, we ignore the region where we
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Left: thermal pressure compared with non-thermal pressure using three different models (black, pink, and green lines, Shaw et al. 2010; Fusco-Femiano & Lapi 2014; Nelson et al. 2014, respectively). Right: measured gas fraction profile azimuthally averaged (red line) and ignoring the merging region (black line), and corrected accounting for the contribution of a non-thermal pressure component enabling to match the cosmic gas fraction at R<sub>200</sub> and R<sub>500</sub>. The horizontal line represents the universal baryon fraction (Planck Collaboration XIII 2016), the vertical lines represents the position of R<sub>500</sub> and R<sub>200</sub>, and the yellow points are the universal baryon fraction depleted by the thermalized gas and by the star fraction. The pink, green, and blue lines represent the gas fraction we get by using different functional forms in order to reduce the observed gas mass fraction to the universal one.</figcaption></figure>

<table><caption>**Table 6.** Comparison between the mass reconstruction at R<sub>200</sub> using the whole surface brightness image and ignoring the merging component.</caption><thead><tr><th>Region</th><th><i>M</i><sub>200</sub> (10<sup>14</sup> <i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>200</sub> (kpc)</th><th><i>M</i><sub>gas,200</sub> (10<sup>14</sup> <i>M</i><sub>☉</sub>)</th><th><i>f</i><sub>gas,200</sub></th></tr></thead><tbody><tr><td>Azimuthal average</td><td>10.7 ± 0.5</td><td>2077 ± 33</td><td>2.54 ± 0.05</td><td>0.237 ± 0.012</td></tr><tr><td>Ignoring the merging region</td><td>10.7 ± 0.3</td><td>2075 ± 17</td><td>2.22 ± 0.02</td><td>0.207 ± 0.006</td></tr></tbody></table>

**Notes.** The columns show: the hydrostatic mass by solving HEE (see Eq. (6)), R<sub>200</sub>, the gas mass obtained by integrating the gas density profile (Eq. (5)), and the gas mass fraction defined by $f_{\text{gas}} = M_{\text{gas}}/M_{\text{tot}}$.

measure this excess in the surface brightness (see red sector in Fig. 13), and repeat our analysis. We show the comparison between the results obtained before and after masking the merging region in Table 6. The hydrostatic mass remains unchanged, but the gas mass decreases, implying that the gas fraction lowers by 17% at R<sub>200</sub>, but it is still greater than the cosmological gas fraction predicted from numerical simulations at these radii. We note that the reconstructed gas fraction is already corrected for the resolved gas clumping using the median density profile; therefore, clumpiness cannot be responsible for the excess gas fraction (Simionescu et al. 2011).

One possibility to explain this overestimate in the gas fraction is the presence of a substantial non-thermal pressure component in the HEE which breaks the hydrostatic equilibrium assumption. We modify the HEE in Eq. (6), by adding an extra pressure component, which we define as non-thermal pressure, and justify as being generated, for example, by unresolved gas turbulence, bulk motion, magnetic field, or asphericity. This non-thermal component can be modelled, in a first approximation, as a constant fraction of the thermal component (Loeb & Mao 1994; Zappacosta et al. 2006). We add this non-thermal pressure term (indicated with the subscript “NT”) in the HEE as $P_{\text{NT}}(r) = \alpha(r)P_{\text{T}}(r)$, where the thermal component has the subscript “T”, and $\alpha(r)$ is a function of radius. The HEE is then modified as

$$ \frac{1}{\rho_g} \left( \frac{dP_T}{dr} + \frac{dP_{NT}}{dr} \right) = -\frac{G}{r^2} (M_T + M_{NT}). \quad (16) $$

By solving the derivatives and readjusting the terms in the equation, we can then write how this propagates into the estimate of the gas mass fraction

$$ \begin{aligned} f_g &= \frac{M_g}{M_T + M_{NT}} = \frac{M_g}{M_T \left(1 + \frac{M_{NT}}{M_T}\right)} \\ &= \frac{f_{g,T}}{1 + \alpha(r) - \frac{P_T r^2}{G M_T \mu m_p n_c} \frac{d\alpha}{dr}} \equiv \beta f_{g,T}, \end{aligned} \quad (17) $$

where $\beta$ is defined as the ratio between the true gas fraction and the measured thermal gas fraction. This means that in the case of $\alpha = \text{const.}$, the real gas fraction is reduced by a factor $1 + \alpha$.

By imposing that the observed cluster gas fraction should match the cosmic value in Eq. (15), and assuming a constant $\alpha$, we require $\alpha = 0.64$ (0.32) at R<sub>200</sub> (R<sub>500</sub>), implying that about 39% (24%) of the total pressure is in the form of a non-thermal component.

In general, $\alpha$ is expected to have a radial dependence. Numerical simulations (e.g., Shaw et al. 2010; Fusco-Femiano & Lapi 2014; Nelson et al. 2014) predict some functional forms for $P_{\text{NT}}/P_T$. We can constrain the parameters of these models, if we consider the radial dependence of $\alpha$ in HEE, by requiring that we should be able to reproduce the expected gas mass fraction at R<sub>500</sub> and R<sub>200</sub>. The errors on the parameters are calculated using Monte Carlo simulations propagating the errors on the gas mass fraction profile, on the measure of R<sub>200</sub>, and on the predicted gas mass fraction points. The non-thermal pressure profiles and the corresponding gas fraction profiles obtained using the above mentioned models are shown in Fig. 16, and in Table 7 we provide the three functional forms adopted and the best-fitting parameters. We observe that already above 200–300 kpc, the non-thermal pressure support plays a very important role in flattening the gas mass fraction profile.
<table><caption>**Table 7.** Model, functional form, and best-fitting parameters for the three models which describe the ratio between non-thermal and thermal pressure support.</caption><thead><tr><th>Model</th><th>Functional form for α = P<sub>NT</sub>/P<sub>T</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>Nelson et al. (2014)</td><td>[a(1 + exp(-((R/R<sub>200</sub>b)<sup>c</sup>)))]<sup>-1</sup> - 1</td><td>0.52 ± 0.02</td><td>0.52 (fix)</td><td>1.23 ± 0.27</td></tr><tr><td>Fusco-Femiano & Lapi (2014)</td><td>a exp(-((1-R/(2R<sub>500</sub>))<sup>2</sup>))</td><td>0.91 ± 0.18</td><td>0.706 ± 0.09</td><td>—</td></tr><tr><td>Shaw et al. (2010)</td><td>a(R/R<sub>500</sub>)<sup>b</sup></td><td>0.63 ± 0.05</td><td>1.17 ± 0.36</td><td>—</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Rescaled pressure (top left) and density (top right) profiles considering the azimuthally averaged, ignoring the merger, and ignoring the merging region and considering the $M_{200,tot}$ and $R_{200,tot}$ required to recover the cosmological gas fraction at the virial radius. We compare these profiles with the Planck envelope (Planck Collaboration Int. V 2013) for pressure, and with the universal density profile (Eckert et al. 2012) for density. Bottom: rescaled entropy and rescaled entropy corrected by the gas mass fraction, before and after correcting for the true total mass.</figcaption></figure>

Finally, by imposing that the total cluster mass $M_{tot}$ is provided from $M_T + M_{NT}$, we can estimate the amount of the hydrostatic bias factor $\beta$ as

$$ \beta = \frac{M_T}{M_{tot}} \Rightarrow M_{tot} = \frac{M_T}{\beta}. \qquad (18) $$

Applying Eqs. (16)–(18), the cosmological gas fraction at $R_{500}$ and $R_{200}$ is obtained by requiring

$$ M_{500,tot} = 10.2 \pm 0.4^{\text{stat.}} \pm 0.4^{\text{syst.}} \times 10^{14} M_{\odot}, $$

$$ M_{200,tot} = 17.3 \pm 0.9^{\text{stat.}} \pm 1.2^{\text{syst.}} \times 10^{14} M_{\odot}. $$

Using this mass estimate corrected by clumpiness and by hydrostatic bias, and the value acquired from the Planck catalogue (Planck Collaboration XXVII 2016) and based on scaling relations, $M_{Y_{SZ},500} = 8.74(\pm0.12) \times 10^{14} M_{\odot}$, we infer a Planck bias of $1 - b = M_{Y_{SZ},500}/M_{500,tot} \approx 0.86$.

## 6.2. Effects of the hydrostatic bias on the rescaled profiles

The correction on the mass propagates to the rescaled profiles, both directly since $R_{500}$ increases shrinking the x-axis, and indirectly since pressure and entropy, as described from Eqs. (8) and (10) respectively, follow a rescaling which is mass dependent.

In Fig. 17, we show the net effect on the thermodynamic rescaled profiles, that can be summarized by the following statements:

- the gas pressure profile is now in agreement with the universal pressure profile (Arnaud et al. 2010) and with the Planck envelope (Planck Collaboration Int. V 2013);

- the gas density profile becomes compatible with the stacked density profile presented in Eckert et al. (2012);

- the gas entropy profile shows the least modification before and after this analysis; the profile becomes slightly steeper, but it is still flat in the outskirts, in agreement with
the expected impact of any non-thermal pressure support (Walker et al. 2012).

Pratt et al. (2010) have shown that in order to reconcile entropy profile with predictions from non-radiative simulations (Voit et al. 2005), the profile has to be corrected by the gas mass fraction $K \Rightarrow K \cdot (E(z)f_{\text{gas}}/f_b)^{2/3}$. Introducing this correction in each entropy profile that we consider (i.e. the azimuthally average profile, the profile ignoring the merging region, and the profile required to recover the cosmological gas fraction at $R_{200}$), we obtain the results shown in Fig. 17. We observe that only when we include the contribution by the non-thermal pressure we obtain a corrected entropy profile that deviates from the numerical predictions, with a flattening above $0.3R_{500}$ suggesting that turbulence, or non-thermal energy at large, has not been yet converted efficiently in heat energy, not allowing the specific entropy of the ICM to rise to the value expected in systems simulated in the absence of non-gravitational processes (e.g., Voit et al. 2005).

## 7. Summary and conclusions

The very accurate background modelling of the XMM-Newton exposures, and the large extension of the SZ signal resolved with Planck have allowed us to combine X-ray and SZ data to study the thermodynamic properties of Abell 2319 over the virial region around $R_{200}$. Moreover, since the data quality is very high, we are able to study the properties of this cluster reaching the virial radius in eight different sectors. This enables us to study the azimuthal variance of the thermodynamic properties of the ICM in this merging system for the first time.

The measured clumpiness shows the presence of the merging component with an increase in its value at intermediate radii (~500 kpc). This excess disappears when we remove the merging regions from the analysis. On the other hand, in the outskirts the clumpiness measured is compatible with the estimated residual clumpiness (Roncarelli et al. 2013). This means that this cluster has no significative infalling clumps at the virial radius.

The gas density profile corrected for the resolved clumpiness is then used to recover other fundamental quantities (Eckert et al. 2015), together with the gas temperature profile that we measure, from the X-ray spectroscopic analysis, with a median relative statistical uncertainty of 2% and with a systematic error that we carefully estimate to be of the order of 4% (median value), and above 15% in the outermost radial bin alone. The exquisite quality of these complementary X-ray and SZ datasets, extending across $R_{200}$, enable us to constrain a NFW hydrostatic mass profile at very high precision ($M_{200} = 10.7 \pm 0.5 \text{stat.} \pm 0.9 \text{syst.} \times 10^{14} M_{\odot}$), achieving a level where systematic errors dominate over the statistical ones.

Due to the merging state of this cluster, the recovered entropy profile is flatter than the one predicted by non-radiative simulations (Voit et al. 2005). We observe the most deviations in the first and last few points: in the centre this is caused by the fact that this cluster is a non-cool core cluster (Cavagnolo et al. 2009) with a flat entropy core of ~75 keV cm$^2$, while some residual non-thermal energy flattens the entropy in the outskirts (Walker et al. 2012).

The pressure profile recovered from SZ data is flatter, and above the $1\sigma$ envelope, than the universal profile measured for an ensemble of objects resolved with Planck (Planck Collaboration Int. V 2013).

The measured gas fraction, corrected by the gas clumpiness using the median density profile, is above the value

predicted from state-of-art hydrodynamical simulations for the preferred cosmological background (Planck Collaboration XIII 2016; Panelles et al. 2013; Gonzalez et al. 2013). Analysing the azimuthal variation of the $f_{\text{gas}}$ profile (see Fig. 14), we observe that it is above the average value only in the sectors most affected by the merger (sectors 1, 2, and 3). When the region with the ongoing merger and with an estimated higher gas mass is excluded from the analysis, the gas fraction drops, but it is still higher than the expectations, indicating a non-negligible contribution from a non-thermal pressure support that we quantify as approximately 39% and 24% of the total pressure at $R_{200}$ and $R_{500}$, respectively.

Once the correction induced by the non-thermal pressure support is propagated through the measurements of $R_{500}$, $K_{500}$, and $P_{500}$, we show that (i) the pressure profile matches the mean behaviour of objects resolved with Planck; (ii) the gas density profile becomes consistent with the stacked profile obtained from Rosat/PSPC observations in Eckert et al. (2012); and (iii) on the contrary, the entropy undergoes a very small change, remaining flatter than the predicted profile.

In forthcoming works, the detailed analysis presented here for A2319 will be extended to the whole X-COP sample (Eckert et al. 2017), providing the first ensembled properties of the ICM at $R_{200}$ and above.

*Acknowledgements.* This research has received funding from the European Union's Horizon 2020 Programme under AHEAD project (grant agreement no. 654215). S.E. acknowledges the financial support from contracts ASI-INAFF I/009/10/0, NARO15 ASI-INAFF I/037/12/0, and ASI 2015-046-R.O.

## References

Ameglio, S., Borgani, S., Pierpaoli, E., & Dolag, K. 2007, MNRAS, 382, 397

Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, ASP Conf. Ser., 101, 17

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012, ApJ, 758, 75

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton University Press)

Bourdin, H., Mazzotta, P., Kozmanyan, A., Jones, C., & Vikhlinin, A. 2017, ApJ, 843, 72

Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, ApJS, 182, 12

Croston, J. H., Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2006, A&A, 459, 1007

De Luca, A., & Molendi, S. 2004, A&A, 419, 837

Diehl, S., & Statler, T. S. 2006, MNRAS, 368, 497

Dolag, K., Bartelmann, M., & Lesch, H. 1999, A&A, 348, 351

Eckert, D., Vazza, F., Ettori, S., et al. 2012, A&A, 541, A57

Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198

Eckert, D., Ettori, S., Coupon, J., et al. 2016, A&A, 592, A12

Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astron. Nachr., 338, 293

Ettori, S., & Molendi, S. 2011, Mem. Soc. Astron. It. Suppl., 17, 47

Ettori, S., De Grandi, S., & Molendi, S. 2002, A&A, 391, 841

Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61

Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68

Ettori, S., Donnarumma, A., Pointecouteau, E., et al. 2013, Space Sci. Rev., 177, 119

Ettori, S., Ghirardini, V., Eckert, D., Dubath, F., & Pointecouteau, E. 2017, MNRAS, 470, L29

Farnsworth, D., Rudnick, L., Brown, S., & Brunetti, G. 2013, ApJ, 779, 189

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, Proc. SPIE, 6270, 62701V

Fusco-Femiano, R., & Lapi, A. 2014, ApJ, 783, 76

Ghizzardi, S., Rossetti, M., & Molendi, S. 2010, A&A, 516, A32

Gonzalez, A. H., Sivanandam, S., Zabludoff, A. I., & Zaritsky, D. 2013, ApJ, 778, 14

Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118

Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775

Kriss, G. A., Cioffi, D. F., & Canizares, C. R. 1983, ApJ, 272, 439

Kuntz, K. D., & Snowden, S. L. 2008, A&A, 478, 575
Lamarre, J.-M., Puget, J.-L., Ade, P. A. R., et al. 2010, A&A, 520, A9

Leccardi, A., & Molendi, S. 2008, A&A, 486, 359

Liu, W., Chiao, M., Collier, M. R., et al. 2017, ApJ, 834, 33

Loeb, A., & Mao, S. 1994, ApJ, 435, L109

Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10

McCammon, D., Almy, R., Apodaca, E., et al. 2002, ApJ, 576, 188

Molendi, S., De Grandi, S., Fusco-Femiano, R., et al. 1999, ApJ, 525, L73

Morandi, A., Ettori, S., & Moscardini, L. 2007, MNRAS, 379, 518

Nagai, D., & Lau, E. T. 2011, ApJ, 731, L10

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25

Oegerle, W. R., Hill, J. M., & Fitchett, M. J. 1995, AJ, 110, 32

Pfrommer, C., Springel, V., Jubelgas, M., & Ensslin, T. A. 2007, in Cosmic Frontiers, eds. N. Metcalfe, & T. Shanks, ASP Conf. Ser., 379, 221

Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1

Planck Collaboration Int. V. 2013, A&A, 550, A131

Planck Collaboration XXIX. 2014, A&A, 571, A29

Planck Collaboration I. 2016, A&A, 594, A1

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration XXVII. 2016, A&A, 594, A27

Planck HFI Core Team 2011, A&A, 536, A4

Panelles, S., Borgani, S., Dolag, K., et al. 2013, MNRAS, 431, 1487

Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85

Read, A. M., Rosen, S. R., Saxton, R. D., & Ramirez, J. 2011, A&A, 534, A34

Roncarelli, M., Ettori, S., Borgani, S., et al. 2013, MNRAS, 432, 3030

Salvetti, D., Marelli, M., Gastaldello, F., et al. 2017, Exp. Astron., 44, 309

Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30

Shaw, L. D., Nagai, D., Bhattacharya, S., & Lau, E. T. 2010, ApJ, 725, 1452

Simionescu, A., Allen, S. W., Mantz, A., et al. 2011, Science, 331, 1576

Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615

Storm, E., Jeltema, T. E., & Rudnick, L. 2015, MNRAS, 448, 2495

Struble, M. F., & Rood, H. J. 1999, ApJS, 125, 35

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comm. Astrophys. Space Phys., 4, 173

Tauber, J. A., Mandolesi, N., Puget, J.-L., et al. 2010, A&A, 520, A1

Tchernin, C., Eckert, D., Ettori, S., et al. 2016, A&A, 595, A42

Tozzi, P., & Norman, C. 2001, ApJ, 546, 63

Vazza, F., Roncarelli, M., Ettori, S., & Dolag, K. 2011, MNRAS, 413, 2305

Vazza, F., Eckert, D., Simionescu, A., Brüggen, M., & Ettori, S. 2013, MNRAS, 429, 799

Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691

Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MNRAS, 364, 909

Walker, S. A., Fabian, A. C., Sanders, J. S., & George, M. R. 2012, MNRAS, 427, L45

Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O'Brien, P. T. 2013, MNRAS, 431, 394

Zappacosta, L., Buote, D. A., Gastaldello, F., et al. 2006, ApJ, 650, 777

Zhuravleva, I., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 428, 3274
<figure><img src="image_19.png" /><figcaption>Fig. A.1. Stacked EPIC radial profiles of 495 blank-sky pointings, sorted in bins of soft-proton contamination inFOV/outFOV. The black data points show observations with low SP contamination (inFOV/outFOV = 1–1.15), whereas the blue points comprise observations that were severely affected by SP contamination (inFOV/outFOV = 1.6–2.0). The left panel shows the stacked profiles obtained when subtracting only the QPB component, while in the right panel, the SP and QC components have been taken into account following Eq. (A.2).</figcaption></figure>

# Appendix A: Non-X-ray background modelling

We developed and calibrated a novel technique to model and subtract the non-X-ray background (NXB). Our approach builds upon the method devised in Tchernin et al. (2016); however, it can be more reliably applied to observations including a significant source emission above 5 keV. Here we describe the main principles of our method and validate it using a large set of blank-sky XMM-Newton pointings.

## A.1. Model

It has long been known that the NXB of XMM-Newton is split into two main components, the quiescent particle background (QPB) and the soft protons (SP). Recently, Salvetti et al. (2017) has analysed almost the complete XMM-Newton archive and showed the presence of an additional stable, low-intensity component within the field of view (FOV) of the MOS2 instrument, whose origin is yet unknown. As described in Sect. 2, a fraction of the area of the MOS detectors is located outside the FOV of the XMM-Newton telescopes. The outFOV area can be used to estimate the QPB level in each observation by rescaling filter-wheel-closed data to the measured outFOV count rate. The remaining inFOV high-energy count rate can then be decomposed into a variable component (SP) and a quiescent part (QC). We can thus describe the remaining NXB as

$$ \text{inFOV} - \text{outFOV} = \text{SP} + \text{QC}, \quad (A.1) $$

where inFOV and outFOV denote the [7–11.5] keV MOS2 count rates measured in the exposed and unexposed areas of the detector, respectively. We restrict the measurement to the MOS2 detector as two of the MOS1 chips have been lost throughout the mission, and the unexposed area of the pn detector is too small for our needs.

Importantly, the SP component is expected to show a different spatial signature on the detector compared to the QPB. Indeed, soft protons, which are funneled towards the detector through the telescope, are more spatially concentrated than the QPB and follow a vignetting curve SP(r) that is different from the vignetting curve of the photons (Kuntz & Snowden 2008),

where r denotes the distance of each pixel from the aim point. Conversely, given that its origin is currently unclear, the spatial distribution of the QC component is unknown. Here we make the hypothesis that this component is flat over the detector.

## A.2. Blank-sky dataset and modelling

To determine the relative contributions of the SP and QC components, we used a large set of 495 XMM-Newton blank-sky pointings, most of which are from the XXL survey (Pierre et al. 2016). Our dataset comprises more than 5 Ms of data. We processed the data using ESAS in the same way as for the A2319 data (see Sect. 2). We estimated the QPB component in each observation by measuring the outFOV count rate and rescaling filter-wheel-closed data. We also compute the high-energy inFOV and outFOV count rates for each observation. We then measured the radial profiles in the [0.7–1.2] keV band of the blank-sky pointings from the aim point to the outermost edge of the pointing in annuli of 30 arcsec width. The detected sources were masked and the QPB was subtracted from the data. As already shown in Tchernin et al. (2016), this procedure results in radial profiles that are on average not flat, which indicates the need of modelling additional components (SP and QC).

We then describe the radial profiles $S_X(r)$ as the sum of the SP and QC components following their respective spatial distributions,

$$ S_X(r) = C + N_{QC} + N_{SP}(\text{inFOV} - \text{outFOV} - \overline{QC})SP(r), \quad (A.2) $$

where C is the sky background intensity at the relevant location, $N_{QC}$ the intensity of the stable QC component, $N_{SP}$ the normalization of the variable SP, and $\overline{QC} = 0.023$ counts s⁻¹ is the mean high-energy count rate of the QC component (Salvetti et al. 2017). We then perform a joint fit on all the measured profiles and optimize for the values of $N_{QC}$ and $N_{SP}$. We then used the best-fit values of $N_{QC}$ and $N_{SP}$ to create 2D models of these components and subtract them from the data.

In Fig. A.1 we show the stacked radial profiles of the full sample. In the left-hand panel we show the stacked profiles obtained when subtracting the QPB component only, whereas
<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** IRAS map (minimum–maximum values in the region within $R_{200}$ are 4.22, 7.77 MJy sr$^{-1}$). The white external circle represents the location of $R_{200}$, while the small circle represents the location of the centre of the cluster.</figcaption></figure>

in the right-hand panel, the SP and QC components have been modelled using the method described above and subtracted from the data. To investigate the dependence of our results on SP contamination, we grouped the data in bins of increasing SP contamination, which we trace using the inFOV/outFOV ratio (Leccardi & Molendi 2008). Observations that were mildly affected by SP contamination exhibit an inFOV/outFOV ratio close to one, whereas heavily contaminated observations show high values of the inFOV/outFOV ratio. The effect of SP contamination is evident in the left-hand panel of Fig. A.1, where the deviations of the stacked profiles from a straight line progressively increase with increasing SP contamination. Conversely, when applying our SP and QC modelling approach, flat profiles are found in all four bins out to the edge of the FOV, indicating that our model accurately reproduces the various NXB components. The excess scatter compared to a straight line is 5%, which we adopt as our systematic uncertainty in the subtraction of the NXB.

## Appendix B: Results of the spectral fitting

In Table B.1, we show the spectral fit results in the analysis described in Sect. 2, indicating the radial extension of the chosen annuli, the C-statistic, the number of the spectral bins, and the reduced C-statistic. We note that this last quantity is always of the order of 1, implying high goodness in the fit.

Since A2319 is located at low galactic latitude, $b = +13.5°$, the choice to leave $n_H$ free to vary is reinforced from the azimuthal variation over the cluster's region of the dust emission as mapped at 100 µm by the InfraRed Astronomical Satellite (IRAS; see Fig. B.1). The map shows that the sectors 5, 6, and 7 are expected to have higher Galactic absorption. Indeed the $n_H$ in the eight considered sectors varies according to Table B.2, with sector 5, 6, and 7 being ~10% above the other sectors.

<table><caption>**Table B.1.** Statistical results of the fitting in the annular regions.</caption><thead><tr><th>Radii<br>arcmin</th><th>C-stat.<br>-</th><th>PHA bins<br>-</th><th>C-stat. reduced<br>-</th><th>Net cts<br>10<sup>3</sup></th><th>SBR<br>-</th><th>n<sub>H</sub><br>10<sup>22</sup> cm<sup>-2</sup></th></tr></thead><tbody><tr><td>0.00–1.05</td><td>2746</td><td>2603</td><td>1.05</td><td>170</td><td>85</td><td>0.075</td></tr><tr><td>1.05–1.63</td><td>2698</td><td>2591</td><td>1.04</td><td>155</td><td>58</td><td>0.078</td></tr><tr><td>1.63–2.18</td><td>2748</td><td>2552</td><td>1.08</td><td>148</td><td>42</td><td>0.081</td></tr><tr><td>2.18–2.74</td><td>2773</td><td>2575</td><td>1.08</td><td>145</td><td>32</td><td>0.077</td></tr><tr><td>2.74–3.32</td><td>2767</td><td>2484</td><td>1.11</td><td>131</td><td>24</td><td>0.078</td></tr><tr><td>3.32–3.98</td><td>2688</td><td>2573</td><td>1.05</td><td>133</td><td>17</td><td>0.081</td></tr><tr><td>3.98–4.65</td><td>2807</td><td>2582</td><td>1.09</td><td>134</td><td>14</td><td>0.079</td></tr><tr><td>4.65–5.37</td><td>2912</td><td>4005</td><td>0.73</td><td>131</td><td>11</td><td>0.075</td></tr><tr><td>5.37–6.14</td><td>2666</td><td>2387</td><td>1.13</td><td>112</td><td>8.6</td><td>0.074</td></tr><tr><td>6.14–6.95</td><td>2811</td><td>2481</td><td>1.13</td><td>101</td><td>6.8</td><td>0.077</td></tr><tr><td>6.95–7.83</td><td>3157</td><td>4949</td><td>0.64</td><td>92</td><td>5.1</td><td>0.074</td></tr><tr><td>7.83–8.85</td><td>3305</td><td>3866</td><td>0.85</td><td>89</td><td>3.6</td><td>0.073</td></tr><tr><td>8.85–10.05</td><td>3697</td><td>6052</td><td>0.61</td><td>82</td><td>2.5</td><td>0.074</td></tr><tr><td>10.05–11.51</td><td>4514</td><td>3868</td><td>1.17</td><td>80</td><td>1.7</td><td>0.076</td></tr><tr><td>11.51–13.10</td><td>4870</td><td>3583</td><td>1.36</td><td>62</td><td>1.2</td><td>0.079</td></tr><tr><td>13.10–15.18</td><td>4893</td><td>3494</td><td>1.40</td><td>46</td><td>0.9</td><td>0.077</td></tr><tr><td>15.18–17.70</td><td>2808</td><td>1844</td><td>1.52</td><td>20</td><td>0.9</td><td>0.121</td></tr><tr><td>17.70–20.63</td><td>2632</td><td>2175</td><td>1.21</td><td>19</td><td>0.6</td><td>0.101</td></tr><tr><td>20.63–24.08</td><td>2098</td><td>1916</td><td>1.09</td><td>12</td><td>0.4</td><td>0.113</td></tr></tbody></table>

**Notes.** The table lists: radial extension, C-statistic, number of spectral bins, reduced C-statistic indicated, net number of photons in the energy band [0.5–11.3] keV, signal-to-background ratio, and best-fit $n_H$.

<table><caption>**Table B.2.** Best-fit $n_H$ in the eight sectors considered.</caption><thead><tr><th>Sector</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td>$n_H$ [10<sup>20</sup> cm<sup>-2</sup>]</td><td>7.65</td><td>7.20</td><td>7.62</td><td>7.87</td><td>8.39</td><td>8.57</td><td>8.41</td><td>7.99</td></tr></tbody></table>

## Appendix C: Comparison with Chandra data

We have analysed two archival *Chandra* observations of the inner region of A2319 (OBSID 15187, with a cleaned exposure time of 75 ks, and OBSID 3231, with 15 ks). We have processed the two *Chandra* ACIS-I observations of A2319 with a standard pipeline based on CIAO 4.9 (Fruscione et al. 2006) and CALDB 4.7.4 to create a new events-2 file which includes filtering for grade, status, bad pixels, and time intervals for anomalous background levels. The background is estimated through blank sky observations. We have extracted the spectra in the same annular regions as for *XMM-Newton*, and fit them in the identical way, leaving the galactic column density $n_H$ free to vary within the range $7–13 \times 10^{20} \text{ cm}^{-2}$. The temperature profiles are compared in Fig. C.1. We observe a good agreement among these spectral measurements, despite the claimed and still debated cross-calibration issue between *Chandra* ACIS and *XMM-Newton* EPIC (see e.g., Schellenberger et al. 2015), in particular in very hot systems ($T > 5$ keV) as A2163. We suggest that leaving free $n_H$ plays a determinant role in adjusting the relative impact of the soft part of the spectra, where most of the observed systematic tension has been reported. In the present case, *Chandra* prefers systematically higher values of $n_H$ (~1.2–1.3 × 10$^{21} \text{ cm}^{-2}$) than *XMM-Newton* (see Table B.1) in all the radial bins. These higher values agree more closely with the column density corrected for molecular hydrogen as suggested in Willingale et al. (2013).
<figure><img src="image_21.png" /><figcaption>Fig. C.1. Comparison between the spectral temperature obtained using Chandra and XMM-Newton. There is a clear excess in the temperature measured by Chandra of the order of 2–3 keV up to 7 arcmin.</figcaption></figure>

## Appendix D: Likelihood for the mass reconstruction

We fit our thermodynamic quantities using the MCMC code emcee (Foreman-Mackey et al. 2013), for which we define a likelihood. We included in the fitting procedure an intrinsic scatter, which is added in quadrature on the error of logarithm of pressure such that $log P \sim log P \pm \sigma_{int}$. By assuming a small value for $\sigma_{int}$ we can write

$$ \sigma_{P,\text{int}} \approx \frac{P \cdot \exp(+\sigma_{\text{int}}) - P \cdot \exp(-\sigma_{\text{int}})}{2} = P \cdot \sinh \sigma_{\text{int}} $$

and summed to the covariance matrix as

$$ \Sigma_{\text{tot}} = \begin{bmatrix} \Sigma_{11} & \Sigma_{12} & \Sigma_{13} & \dots & \Sigma_{1n} \\ \Sigma_{21} & \Sigma_{22} & \Sigma_{23} & \dots & \Sigma_{2n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \Sigma_{n1} & \Sigma_{n2} & \Sigma_{n3} & \dots & \Sigma_{nn} \end{bmatrix} + \begin{bmatrix} \sigma_{P_1,\text{int}}^2 & 0 & 0 & \dots & 0 \\ 0 & \sigma_{P_2,\text{int}}^2 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \sigma_{P_n,\text{int}}^2 \end{bmatrix}, $$

where $\Sigma_{i,j}$ is the covariance matrix on the measured Planck pressure profile.

The intrinsic scatter is also propagated to the variance on temperature profile, added in quadrature to the measured errors:

$$ \sigma_{\text{tot}}^2 = \sigma_{\text{T}}^2 + \sigma_{\text{T,int}}^2 $$

with

$$ \sigma_{\text{T,int}} = \frac{P_{\text{model}}}{n_{\text{model}}} \sigma_{\text{P,int}} = T_{\text{model}} \cdot \sigma_{\text{P,int}}. $$

We recall that in general the likelihood is defined as

$$ \mathcal{L} = \frac{1}{\sqrt{2\pi}\sigma^2} \exp(-\chi^2/2), $$

so that

$$ \log \mathcal{L} = -0.5(\chi^2 + \log \sigma^2 + \log(2\pi)), $$

where the last term is a constant, and therefore is usually ignored while maximizing the likelihood, but the term with $\log \sigma^2$ is not. Finally, by using the subscript “m” or “o” to describe model predicted or observed quantities, respectively, we can explicitly write the logarithm of the likelihood we use to fit

$$ \log \mathcal{L} = -0.5 \left[ (P - P_m) \Sigma_{\text{tot}}^{-1} (P - P_m)^T + n \log (\det (\Sigma_{\text{tot}})) \right] \\ - 0.5 \sum_{i=1}^{n} \left[ \frac{(T_i - T_{m,i})^2}{\sigma_{T,i}^2 + \sigma_{T,int}^2} + \log (\sigma_{T,i}^2 + \sigma_{T,int}^2) \right] \\ - 0.5 \left[ \sum_{i=1}^{n} \frac{(\epsilon - \epsilon_{m,i})^2}{\sigma_{\epsilon,i}^2} \right]. $$

We point out that this method is independent of the method used to compute $P_{model}$ and $T_{model}$, meaning that this kind of approach is valid both for the forward and backward methods.

## Appendix E: Thermodynamic quantities in azimuthal sectors

The procedure described in Sects. 2–4 are applied on each azimuthal sector. In summary, we deproject surface brightness into density using the multiscale technique on the mean profile, we deproject the Comptonization parameter to retrieve pressure, and we calculate the temperature in six spectral annuli. We then apply the backward approach to these thermodynamic quantities in order to find the parameters of a NFW mass model which best reproduce the observables. We compare the observed and reconstructed from the best-fit mass model pressure and temperature profiles sector by sector in Figs. E.1 and E.2, respectively. We observe that the only sectors with an evident discrepancy are the ones disturbed the most by the merger event, i.e., sectors 1, 2, and 3.

Similarly to what is done in Sect. 4.5, we compare the entropy profile reconstructed by the NFW backward best fit with the entropy recovered from X-ray spectroscopy ($K = kT/n_e^{2/3}$), and with the entropy recovered by combining X-ray density and SZ pressure ($K = P/n_e^{5/3}$); this sector-by-sector comparison is shown in Fig. E.3.
<figure><img src="image_22.png" /><figcaption>**Fig. E.1.** Comparison of the observed pressure profile with that reconstructed by the NFW *backward* best fit.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. E.2.** Comparison of the observed 2D temperature profile with the one reconstructed by the NFW *backward* best fit.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. E.3.** Comparison between the entropy profile reconstructed by the NFW *backward* best fit with the entropy coming from the combination of X-ray and SZ and just using X-ray spectral results.</figcaption></figure>