© ESO 2018

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).

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 suchas 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 Ω_Λ = 0.7, Ω _m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹, and . At the redshift of A2319, 1 arcmin corresponds to approximatively 64.9 kpc. Uncertainties are provided at the 1σ confidence level.

In the following, we refer to and plot as reference characteristic radii, R₅₀₀ = 1368 kpc and R₂₀₀ = 2077 kpc, that are defined at the overdensities of Δ = 500 and 200, respectively, with respect to the critical value 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 MOS 1, MOS 2, and pn and mosaicked all seven observations into one single image.

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 orderto 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 ± 0.06) × 10⁻⁴ cts s⁻¹ arcmin⁻², 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 ± 0.05 × 10⁻¹⁵ erg s⁻¹ cm⁻² arcmin⁻²).

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 MOS 2 (since the combined image was in units of MOS 2), 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σ, and mean relative deviation of 5% up to the virial radius (see Fig. 3).

Fig. 1 Mosaicked and Voronoi tessellated image of A2319 in the energy band [0.7–1.2] keV correctedfor 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 R200.

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 R500 and R200, respectively.

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 R500 and R200, respectively.

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 withstrong 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. 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 × 10⁻¹⁵ erg s⁻¹ cm⁻² arcmin⁻² 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 tothe covered areas.

We fixed the particle and sky background. The cluster emission was modelled with a thermal component absorbed from our 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²⁰ and 12.8 × 10²⁰ cm⁻², 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 (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): , where 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 . Values of the temperature and metal abundance in each shell are then obtained as , where 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 (2)

where the integral is computed along the line of sight, ℓ, at the radius, r, from the centre; σ_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 Δ θ∕θ₅₀₀ = 0.2. The local background offset was estimated from the area surrounding the cluster beyond 5 × θ₅₀₀ = 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 even at the moderate angular resolution of the Planck survey, and its SZ signal extends well beyond R₅₀₀ 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

where Σ 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.

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 R500 and R200, respectively.

Fig. 5 Comptonization profile extracted from the SZ map. The vertical dotted and dashed line indicates the location of R500 and 2R500, respectively.

Fig. 6 Planck correlation matrix ρX,Y for the unbinned Comptonization parameter profile.

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

(3)

• the gas entropy

(4)

• the gas mass

(5)

where the gas mass density ρ_g = (n_e + n_p)m_uμ with m_u being the atomic mass unit and μ ≈ 0.6 the mean molecular weight in amu;

• the hydrostatic gravitating mass:

(6)

where G is the gravitational constant, and the gas pressure P_g satisfies the ideal gaslaw ρ_gkT∕(μm_u) = P_g. The gas mass fraction is then defined as f_gas = M_gas∕M_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 measures the bias that affects the reconstruction of the gas density from the X-ray emission, which is directly proportional to . Since we areconsidering 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 of about 1.1 at R₂₀₀.

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 ≈ 18.4 kpc, for MOS¹; see also Read et al. (2011)). This implies that clumped structures below this scale might still bias our measured thermodynamic quantities.

Fig. 7 Clumping factor radial profile for both techniques, onion peeling (blue) and multiscale (red). The pink area represents the 1σ 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 R500 and R200, respectively.

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.4and (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 respectto 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 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%.

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σ confidence region around the backward result. The vertical dotted and dashed line marks the position of R500 and R200, respectively.

Fig. 9 Temperature and abundance profile adopting different techniques; this provides an estimate of the systematic error affecting our measurement.

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), (7)

where (8)

and x = R∕R₅₀₀; γ, α, and β are the central slope, the intermediate slope, and the outer slope defined by a scale parameter r_s = R₅₀₀∕c₅₀₀ (R ≪ r_s, R ~ r_s and R ≫ r_s respectively); and P₀ is the normalization. The values of R₅₀₀ and M₅₀₀ 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.

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₅₀₀ 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₂₀₀, 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.

Fig. 10 Rescaled pressure profile in units of R500. The grey dashed line represents the universal pressure profile. The blue and the pink lines representsthe best fits using the functional form introduced by Nagai et al. (2007) done on PNFW and PSZ, respectively. The dotted and dashed vertical line represents the position of R500 and R200, respectively. In the bottom panel we show the ratio of PSZ, PX, and PNFW with the universal pressure profile (Arnaud et al. 2010).

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 ina forthcoming publication), and constrain its two parameters, concentration and scale radius (or R₂₀₀), using both the projected temperatureprofile 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 yieldsdirectly 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₂₀₀, direct spectral estimates of the temperature are limited to regions below R₅₀₀, 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.

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₂₀₀ 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₅₀₀ and R₂₀₀ 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₅₀₀ and R₂₀₀, respectively, implying that the reference values for the hydrostatic mass are, at R₅₀₀ and R₂₀₀, respectively:

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 emissivityprofile 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 R500 and R200, respectively.

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) (9)

The central entropy (k₀) 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₀ < 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) (10)

where f_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: (11)

In orderto accommodate the flattening of the entropy profile observed in many disturbed systems, we add a constant to a simple power law: (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).

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 KNFW data using Eq. (9). Bottom: entropy profiles rescaled by K500. The dashed magenta line represents the best fit obtained on the KNFW data using Eq. (12). The green lines represents 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 R500 and R200, respectively.

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.

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 6shows the cold front that was detected in Ghizzardi et al. (2010) and located in the SE region, about 200 kpc ≈ 3 arcmin from thecluster 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₂₀₀ (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.

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 R200. The red sector represents the region which shows a clear excess in the residual map.

5.2 Azimuthal scatter and clumpiness

The azimuthal scatter of the recovered thermodynamic quantities is defined at each radius r as (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 σ_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₅₀₀ and R₂₀₀ where the azimuthal scatter reaches a minimum. This point suggests the radial 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 (14)

where r = R∕R₂₀₀; σ is the azimuthal scatter of the density n or of the Comptonization parameter y; and σ₀ and r₀ are estimated from simulations (Roncarelli et al. 2013)

• (σ₀, r₀) = (16.02, 5.87) for X-ray density,

• (σ₀, r₀) = (2.83, 8.25) 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₂₀₀ decreases to 1.05.

Fig. 14 Topleft: 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 R500 and R200, 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 R500 and R200, respectively.

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 R500 and R200, 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 R500 and R200, 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 R500 and R200, respectively.

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), (15)

where Ω_b and Ω _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 Ω _b = 0.045 and Ω _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₅₀₀ and R₂₀₀, 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₅₀₀ and R₂₀₀, 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₅₀₀ (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 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₂₀₀, 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_NT(r) = α(r)P_T(r), where the thermal component has the subscript “T”, and α(r) is a functionof radius. The HEE is then modified as (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 (17)

where β is defined asthe ratio between the true gas fraction and the measured thermal gas fraction. This means that in the case of α =const., the real gas fraction is reduced by a factor 1 + α.

By imposing that the observed cluster gas fraction should match the cosmic value in Eq. (15), and assuming a constant α, we require α = 0.64 (0.32) at R₂₀₀ (R₅₀₀), implying that about 39% (24%) of the total pressure is in the form of a non-thermal component.

In general, α 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_NT∕P_T. We can constrain the parameters of these models, if we consider the radial dependence of α in HEE, by requiring that we should be able to reproduce the expected gas mass fraction at R₅₀₀ and R₂₀₀. 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₂₀₀, 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.

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 β as (18)

Applying Eqs. (16)–(18), the cosmological gas fraction at R₅₀₀ and R₂₀₀ is obtained by requiring

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, , we infer a Planck bias of .

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 R200 and R500. The horizontal line represents the universal baryon fraction (Planck Collaboration XIII 2016), the vertical lines represents the position of R500 and R200, 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.

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₅₀₀ 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 . 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₂₀₀), 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₅₀₀ 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).

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 M200,t ot and R200,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.

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₂₀₀. 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₂₀₀, enable us to constrain a NFW hydrostatic mass profile at very high precision (M₂₀₀ = 10.7 ± 0.5^stat. ± 0.9^syst. × 10¹⁴ M_⊙), 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², 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σ 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; Planelles et al. 2013; Gonzalez et al. 2013). Analysing the azimuthal variation of the f_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₂₀₀ and R₅₀₀, respectively.

Once the correction induced by the non-thermal pressure support is propagated through the measurements of R₅₀₀, K₅₀₀, and P₅₀₀, 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₂₀₀ 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-INAF I/009/10/0, NARO15 ASI-INAF I/037/12/0, and ASI 2015-046-R.0.

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 MOS 2 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 decomposedinto a variable component (SP) and a quiescent part (QC). We can thus describe the remaining NXB as (A.1)

where inFOV and outFOV denote the [7–11.5] keV MOS 2 count rates measured in the exposed and unexposed areas of the detector, respectively. We restrict the measurement to the MOS 2 detector as two of the MOS 1 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, (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 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 ofN_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 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.

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).

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 toleave 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.

Fig. B.1 IRAS map (minimum–maximum values in the region within R200 are 4.22, 7.77 MJy sr−1). The white external circle represents the location of R200, while the small circle represents the location of the centre of the cluster.

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 × 10²⁰ cm⁻². 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²¹ cm⁻²) 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).

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.

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 includedin the fitting procedure an intrinsic scatter, which is added in quadrature on the error of logarithm of pressure such that log P ~ logP ± σ_int. By assuming a small value for σ_int we can write

and summed to the covariance matrix as

where Σ_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:

with

We recall that in general the likelihood is defined as

so that

where the last term is a constant, and therefore is usually ignored while maximizing the likelihood, but the term with log σ² 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

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.

Fig. E.1 Comparison of the observed pressure profile with that reconstructed by the NFW backward best fit.

Fig. E.2 Comparison of the observed 2D temperature profile with the one reconstructed by the NFW backward best fit.

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.

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 (), and with the entropy recovered by combining X-ray density and SZ pressure (); this sector-by-sector comparison is shown in Fig. E.3.

References

All Tables

All Figures

	Fig. 1 Mosaicked and Voronoi tessellated image of A2319 in the energy band [0.7–1.2] keV correctedfor 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 R200.
In the text

	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 R500 and R200, respectively.
In the text

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 R500 and R200, respectively.

In the text

	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 R500 and R200, respectively.
In the text

	Fig. 5 Comptonization profile extracted from the SZ map. The vertical dotted and dashed line indicates the location of R500 and 2R500, respectively.
In the text

	Fig. 6 Planck correlation matrix ρX,Y for the unbinned Comptonization parameter profile.
In the text

	Fig. 7 Clumping factor radial profile for both techniques, onion peeling (blue) and multiscale (red). The pink area represents the 1σ 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 R500 and R200, respectively.
In the text

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σ confidence region around the backward result. The vertical dotted and dashed line marks the position of R500 and R200, respectively.

In the text

	Fig. 9 Temperature and abundance profile adopting different techniques; this provides an estimate of the systematic error affecting our measurement.
In the text

Fig. 10 Rescaled pressure profile in units of R500. The grey dashed line represents the universal pressure profile. The blue and the pink lines representsthe best fits using the functional form introduced by Nagai et al. (2007) done on PNFW and PSZ, respectively. The dotted and dashed vertical line represents the position of R500 and R200, respectively. In the bottom panel we show the ratio of PSZ, PX, and PNFW with the universal pressure profile (Arnaud et al. 2010).

In the text

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 emissivityprofile 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 R500 and R200, respectively.

In the text

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 KNFW data using Eq. (9). Bottom: entropy profiles rescaled by K500. The dashed magenta line represents the best fit obtained on the KNFW data using Eq. (12). The green lines represents 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 R500 and R200, respectively.

In the text

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 R200. The red sector represents the region which shows a clear excess in the residual map.

In the text

Fig. 14 Topleft: 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 R500 and R200, 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 R500 and R200, respectively.

In the text

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 R500 and R200, 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 R500 and R200, 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 R500 and R200, respectively.

In the text

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 R200 and R500. The horizontal line represents the universal baryon fraction (Planck Collaboration XIII 2016), the vertical lines represents the position of R500 and R200, 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.

In the text

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 M200,t ot and R200,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.

In the text

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).

In the text

	Fig. B.1 IRAS map (minimum–maximum values in the region within R200 are 4.22, 7.77 MJy sr−1). The white external circle represents the location of R200, while the small circle represents the location of the centre of the cluster.
In the text

	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.
In the text

	Fig. E.1 Comparison of the observed pressure profile with that reconstructed by the NFW backward best fit.
In the text

	Fig. E.2 Comparison of the observed 2D temperature profile with the one reconstructed by the NFW backward best fit.
In the text

	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.
In the text