Multi-probe analysis of the galaxy cluster CL J1226.9+3332: Hydrostatic mass and hydrostatic-to-lensing bias

The precise estimation of the mass of galaxy clusters is a major issue for cosmology. Large galaxy cluster surveys rely on scaling laws that relate cluster observables to their masses. From the high resolution observations of ~ 45 galaxy clusters with NIKA2 and XMM-Newton instruments, the NIKA2 SZ Large Program should provide an accurate scaling relation between the thermal Sunyaev-Zel'dovich effect and the hydrostatic mass. In this paper, we present an exhaustive analysis of the hydrostatic mass of the well known galaxy cluster CL J1226.9+3332, the highest-redshift cluster in the NIKA2 SZ Large Program at z = 0.89. We combine the NIKA2 observations with thermal Sunyaev-Zel'dovich data from NIKA, Bolocam and MUSTANG instruments and XMM-Newton X-ray observations and test the impact of the systematic effects on the mass reconstruction. We conclude that slight differences in the shape of the mass profile can be crucial when defining the integrated mass at R500, which demonstrates the importance of the modeling in the mass determination. We prove the robustness of our hydrostatic mass estimates by showing the agreement with all the results found in the literature. Another key information for cosmology is the bias of the masses estimated assuming hydrostatic equilibrium hypothesis. Based on the lensing convergence maps from the Cluster Lensing And Supernova survey with Hubble (CLASH) data, we obtain the lensing mass estimate for CL J1226.9+3332. From this we are able to measure the hydrostatic-to-lensing mass bias for this cluster, that spans from 1 - bHSE/lens ~ 0.7 to 1, presenting the impact of data-sets and mass reconstruction models on the bias.


Introduction
Galaxy clusters are formed by gravitational collapse at the last step of the hierarchical structure formation process (Kravtsov & Borgani 2012). Thus, they are tracers of the large-scale structure formation physics. Their abundance in mass and redshift is sensitive to the initial conditions in the primordial Universe and its expansion history and matter content (Huterer et al. 2015). Therefore, galaxy clusters are probes of the underlying cosmology (Allen et al. 2011).
Large catalogs of clusters (e.g. Planck Collaboration XVIII. 2015; Rykoff et al. 2016;Adami et al. 2018;Bleem et al. 2020) have enabled, in the past few years, to constrain cosmological parameters. Nevertheless, these results show some tension with respect to the cosmology obtained from the Cosmic microwave background (CMB) power spectrum analysis (Planck Collaboration XXI. 2013;Planck Collaboration XXIV. 2015;Salvati et al. 2018), in line with a more general problem in that early-and late-Universe probes are giving different results (Verde et al. 2019). Cluster-based cosmological analyses rely on their distribution in mass and one source of the discrepancy may be the inaccuracy on those mass estimates (e.g. Pratt et al. 2019;Salvati et al. 2020). miren.munoz@lpsc.in2p3.fr About 85% of the total mass of clusters of galaxies is composed of dark matter. The remaining 15% corresponds to the hot intra-cluster medium (ICM) and the galaxies in the cluster. For this reason, most of the mass content of the clusters is not directly observable and it needs to be estimated either from the gravitational potential reconstruction or via scaling relations that link cluster observables to their masses (Pratt et al. 2019). The gravitational potential of clusters can be inferred from their gravitational lensing effect on background sources (Bartelmann 2010), from the dynamics of member galaxies (Biviano & Girardi 2003;Aguado-Barahona et al. 2022) or from the combination of the thermodynamical properties of the gas in the ICM under hydrostatic equilibrium (HSE) hypothesis. Scaling relations that allow to recover the mass have been measured and calibrated for various cluster observables at different wavelengths: in the optical and infrared domains the galaxy member richness (Andreon & Hurn 2010), in X-ray the emitted luminosity of the hot gaseous ICM (Pratt et al. 2009) and in millimeter wavelengths the amplitude of the Sunyaev-Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1972;Planck Collaboration XI. 2011) proportional to the thermal energy in the ICM.
The thermal SZ (tSZ) effect (Sunyaev & Zeldovich 1972), corresponding to the spectral distortion of the CMB photons due to the inverse Compton scattering on the hot thermal electrons of Article number, page 1 of 23 arXiv:2209.07460v1 [astro-ph.CO] 15 Sep 2022 the ICM, is an excellent way of detecting clusters as its amplitude is not affected by cosmological dimming. Large surveys, in particular the Planck satellite observations (Planck Collaboration XVIII. 2015) and the ground-based Atacama Cosmology Telescope (ACT, Hilton et al. 2018;Hilton et al. 2021) and South Pole Telescopes (SPT, Bleem et al. 2020), have obtained large SZ-detected galaxy cluster catalogs. However, they need to rely on the aforementioned SZ-mass scaling relation in order to carry out cosmological analyses. So far, most of the SZ-mass scaling relations have been mainly determined from low-redshift (z < 0.5) cluster samples with masses obtained from X-ray observations Planck Collaboration XI. 2011). Other scaling relations, with optical data for example (Saro et al. 2017), are also used.
In any case, it is essential to study the redshift evolution of these scaling relations, as they would impact the cosmological results (Salvati et al. 2020). When building scaling relations, another key aspect is the impact of the assumptions on which the mass computations rely, such as the HSE hypothesis. High spatial resolution cluster observations can assess these assumptions by studying the impact of the dynamical states of clusters on the mass-observable scaling relations.
The NIKA2 SZ Large program Perotto et al. 2022), described in Sect. 3.1.1, seeks to address the aforementioned issues. The work presented in this paper constitutes the third analysis of a cluster in the NIKA2 SZ Large Program. The first analysis on PSZ2 G144.83+25.11 comprised a science verification study, as well as the proof of the impact of substructures in the reconstruction of the physical cluster properties . The second, the worst-case scenario for the NIKA2 SZ Large Program, analysed the ACT-CL J0215.4+0030 galaxy cluster, proving the quality of NIKA2 camera Bourrion et al. 2016;Calvo et al. 2016;Perotto et al. 2020) in the most challenging case of a high-redshift and low mass cluster (Kéruzoré et al. 2020). In this work we present a study on the HSE mass of the CL J1226.9+3332 galaxy cluster, the impact of different systematic effects on the recovered mass, as well as a thorough comparison to previous works. Moreover, we re-estimate its lensing mass and compare it to the mass obtained under the hydrostatic assumption, therefore allowing us to compute the hydrostatic-to-lensing mass bias. This paper is organized as follows. In Sect. 2 we present a summary of the results in the literature for CL J1226.9+3332. In Sect. 3 we describe the observations of the ICM with NIKA2 and XMM-Newton instruments. The reconstruction of the thermodynamical properties of the ICM is presented in Sect. 4. We detail the pressure profile reconstruction from NIKA2 maps, accounting for systematic effects due to the data reduction process and point source contamination. We compare these profiles to previous results as well as to the X-ray pressure profile. In Sect. 5 we describe the HSE mass estimates from the combination of SZ and X-ray data and test the robustness of the results against data processing and modeling effects. In Sect. 6 we present the lensing masses that we estimate for CL J1226.9+3332. All the results are put together in Sect. 7 for discussion, where we compute the hydrostatic-to-lensing mass bias for this cluster. The conclusions are given in Sect. 8.
Throughout this paper we assume a flat ΛCDM cosmology with H 0 = 70 km/s/Mpc and Ω m,0 = 0.3. With this assumption, 1 arcmin corresponds to a distance of 466 kpc at the cluster redshift, z = 0.89.

Previous observations
Since the first SZ observations with BIMA (Joy et al. 2001), the projected morphology of CL J1226.9+3332 appeared to have a quite circular symmetry. Nevertheless, the combination of XMM-Newton and Chandra X-ray data (Maughan et al. 2007) showed a region, at ∼ 40" to the south-west of the X-ray peak, with much higher temperature than the average in the ICM. This substructure was also confirmed by posterior SZ analyses with MUSTANG (Korngut et al. 2011) and NIKA (Monfardini et al. 2011;Adam et al. 2015;. Romero et al. (2018), hereafter R18, performed a study combining SZ data from NIKA, MUSTANG and Bolocam instruments. Their different capabilities enabled to probe different angular scales in the reconstruction of ICM properties and agreed with a nonrelaxed cluster core description for CL J1226.9+3332. In this work we will make use of the pressure profiles obtained from NIKA, MUSTANG and Bolocam data summarized in Table 2 in Romero et al. (2018).
Lensing data from the Cluster Lensing And Supernova survey with Hubble (CLASH, Zitrin et al. 2015), as well as the galaxy distribution in the cluster (Jee & Tyson 2009), agree on the existence of a main clump centered on the BCG and a secondary clump on its south-west.
However, this second region does not appear as a structure in X-ray surface brightness (Maughan et al. 2007). One hypothesis presented in Jee & Tyson (2009) suggests that the mass of the south-western galaxy group is not big enough to be observed as an X-ray overdensity. Motivated by the slight elongation of the X-ray peak toward the southwest, Jee & Tyson (2009) also hypothesize that the two-halo system is being observed after the less massive cluster has passed through the central one. A previous study (Maughan et al. 2004) showed also a region of cooler emission in the west side of the BCG, that is, in the north of the mentioned hot region. This was seen using Chandra data and it was explained as a possible infall of some cooler body. Additionally, from the diffuse radio emission analysis with LOFAR data, Di Gennaro et al. (2020) showed that CL J1226.9+3332 hosts the most distant radio halo discovered to date: a radio emission with a size of 0.7 Mpc that follows the thermal gas distribution. In brief, CL J1226.9+3332 shows evidence of disturbance in the core, but a relaxed morphology at large scales. Full grey diamonds represent HSE masses from the combination of SZ and X-ray data and empty ones correspond to X-ray-only results. Pink squares show the SZ-only mass assuming Virial relation, purple circles are dynamical mass estimates and brown stars correspond to lensing M 500 .

The mass of CL J1226.9+3332
Regarding the mass of CL J1226.9+3332, which constitutes the main topic in this study, we present here the results obtained in previous works (summarized in Table A.1). These masses have not been homogenised, nor scaled to the same cosmology and are the values extracted directly from different analyses. We define M ∆ as the mass of the cluster inside a sphere of radius R ∆ . R ∆ is the radius at which the mean mass density of the cluster is ∆ times the critical density of the Universe at its redshift, ρ c = 3H(z) 2 /(8πG) where H(z) is the Hubble function.
We display in Fig. 1 the different M 500 estimates found in the literature. Grey diamonds with error bars correspond to the HSE mass estimates. We distinguish the HSE masses obtained from the combination of SZ and X-ray data (full diamonds) and the X-ray-only results (empty diamonds). The mass given by Planck Collaboration XVIII. (2015) is considered here a SZ+X result, but it is important to keep in mind that this mass was ob-tained applying a scaling relation (derived from X-ray data) to the Planck SZ signal. The empty pink squares show the M 500 assuming the Virial relation and using only SZ data (Mroczkowski 2011(Mroczkowski , 2012. The purple circles are the dynamical masses from Aguado- Barahona et al. (2022) and the brown stars the lensing estimates from Merten et al. (2015) and Sereno (2015). We decide not to present in the same figure the projected masses, as it would be misleading to compare them to the masses integrated in a sphere. The figure shows that both HSE and lensing masses among them vary more than 40% from one analysis to another.
All these mass estimates for CL J1226.9+3332 are hindered by systematic effects, which are difficult to deal with. Comparing properly masses obtained from different observables, methods or modeling approaches is crucial, but very challenging. Moreover, as the shape of the HSE mass profile varies depending on the data and the analysis procedure that is considered, the value of R 500 is not the same for all estimates presented in Fig. 1. Comparisons are thus delicate due to the correlation between the mass and the radius at which it is estimated. As mentioned, an accurate knowledge of the mass of galaxy clusters will be essential for cosmological purposes (Pratt et al. 2019). This motivates the study in this paper, which follows previous work in Ferragamo et al. (2022).

NIKA2 and the SZ Large Program
The NIKA2 camera Bourrion et al. 2016;Calvo et al. 2016) is a millimeter camera operating at the IRAM 30-meter telescope on Pico Veleta (Sierra Nevada, Spain). It observes simultaneously in two bands centered at 150 and 260 GHz with three arrays of kinetic inductance detectors (KIDs, Day et al. 2003;Shu et al. 2018), two operating at 260 GHz and one at 150 GHz. These frequency bands are adapted to detect the tSZ effect characteristic distortion, which shows a decrement of the CMB brightness at frequencies lower than ∼ 217 GHz and an increment at higher frequencies. As described in the instrument performance analysis review (Perotto et al. 2020), NIKA2 maps the sky with a 6.5' field of view and high angular resolution, 17.6" and 11.1" full width at half maximum at 150 and 260 GHz, respectively. These capabilities, combined with its high sensitivity (9 mJy·s 1/2 at 150 GHz and 30 mJy·s 1/2 at 260 GHz), enable to study the ICM of galaxy clusters in detail, as it has already been proved in previous works Kéruzoré et al. 2020). Substructures in clusters as well as contaminant sources can also be identified.
Part of the NIKA2 guaranteed time was allocated to the NIKA2 SZ Large Program. This program is a high angular resolution follow-up of ∼ 45 galaxy clusters detected with Planck and ACT (Planck Collaboration XVIII. 2015; Hasselfield et al. 2013). These clusters have been chosen at intermediate to high redshift ranges (0.5 ≤ z ≤ 0.9) and covering a wide range of masses (3 ≤ M 500 /10 14 M ≤ 10) estimated according to their tSZ signal and mass-observable scaling relations. The improvement in resolution with respect to the previous instruments, going below the arcminute scale, allows us to resolve distant clusters for which the apparent angular sizes are small. The SZ Large Program also benefits from the X-ray observations obtained with the XMM-Newton satellite. From the combination of tSZ and X-ray data, the NIKA2 SZ Large Program will be able to reestimate precise HSE masses of the galaxy clusters in the sample and improve the current SZ-mass scaling relations. A more detailed explanation of the NIKA2 SZ Large Program is given in Mayet et al. (2020); Perotto et al. (2022).

NIKA2 observations and maps of CL J1226.9+3332
CL J1226.9+3332 was observed for 3.6 hours during the 15th NIKA2 science-purpose observation campaign (13-20 February 2018) as part of the NIKA2 Guaranteed Time (under project number 199-16) at the IRAM 30-m telescope. The data consists of 36 raster scans of 8 × 4 arcminutes in series of four scans with angles of 0, 45, 90, and 135 degrees with respect to the rightascension axis. The scans were centered at the XMM-Newton Xray peak, (R.A., Dec.) J2000 = (12h26m58.08s, +33d32m46.6s). The mean elevation of the scans is of 58.51 • . Data at 150 and 260 GHz were acquired simultaneously.
The raw data were calibrated and reduced following the baseline method developed for the performance assessment of NIKA2 (Perotto et al. 2020) and extended to diffuse emission as in Kéruzoré et al. (2020) and Rigby et al. (2021). The data for the two frequency bands are processed independently. Common modes of most-correlated detectors are estimated for data outside a circular mask covering the cluster. These modes are subtracted from the raw data and a new mask is defined at the positions where the significance of the signal is above a threshold, S/N > 3 in this case. At the next iteration the common modes are estimated outside the new mask and we repeat the procedure until convergence. We chose as initial mask a 2 arcminute diameter disk centered at the center of the scans.
We present in Fig. 2 the resulting NIKA2 surface brightness maps at 150 and 260 GHz for CL J1226.9+3332. Black contours indicate significance levels starting from 3σ with a 3σ spacing. The map at 150 GHz (left panel) shows the cluster as a negative decrement with respect to the background. This is the characteristic signature of the tSZ effect at frequencies lower than 217 GHz. In this map we also identify positive sources that can compensate the negative tSZ signal of the cluster, as it is the case for the central southeastern source. Moreover, in the 150 GHz map we observe an elongation of the tSZ peak towards the south west. Similar structures were found by Maughan et al. (2007); Korngut et al. (2011);Adam et al. (2015); Zitrin et al. (2015) and Jee & Tyson (2009) as mentioned in Sect. 2. The 260 GHz map, in the right panel, is dominated by the signal of the point sources and can be used to identify them. Nonetheless, tSZ signal of the cluster is not detected at this frequency. This is expected, since the integrated tSZ signal is about three times weaker in the NIKA2 band centered at 260 GHz than in the 150 GHz one.

Estimation of noise residuals in the NIKA2 maps
The residual noise in the final 150 GHz NIKA2 map of CL J1226.9+3332 needs to be quantified for the reconstruction of the pressure profile of the cluster. It is usually estimated on null maps, also known as jackknives (JK), by computing halfdifferences of two statistically equivalent sets of scans in order to eliminate the astrophysical signal and recover the residuals. We have developed two different noise estimates to evaluate possible systematic bias and uncertainties. The so called "angle order" (AO) noise map is computed from the half-differences of scans observed with the same angle with respect to the rightascension axis. This ensures that signal residuals from differential filtering along the scan direction are minimized in the null maps. Alternatively, the so-called "time order" (TO) noise map is calculated from the half-differences of consecutive scans. This  . Power spectra of noise map estimates for the NIKA2 150 GHz data: in pink and black the spectra for the JK maps estimated with angle ordered and time ordered scans, respectively. In blue and green we show the power spectra of the different residual maps for the best-fit models shown in Fig. 5. The grey regions represent the NIKA2 instrumetal limits given by the field of view (for small angular frequencies) and the beam FWHM (for large angular frequencies). minimizes the time dependent effects that may be induced by atmospheric residual fluctuations. We present in Fig. 3 in pink and black the power spectra of the AO and TO null maps at 150 GHz for CL J1226.9+3332, respectively. The AO null map has a flatter power spectrum for small wave numbers, meaning that it contains less large-scale correlated noise than the TO null map. This would suggest that the TO null map might be affected by signal or atmospheric residuals differently filtered for each scanning angle. The power spectra shown in Fig. 3 have been used to compute 1000 Monte Carlo noise realizations to estimate the pixel-pixel noise covariance matrices used in Sect. 4.1.1 (following the method developed in Adam et al. 2016).

Pipeline Transfer function estimation
The filtering induced by the observations and the data reduction process on the cluster signal needs to be evaluated to be accounted for when estimating the pressure profile of the cluster. The transfer function is calculated by repeating the data processing steps discussed above on a simulation of the cluster signal. It is computed in Fourier space as the ratio of the power spectrum of the simulation filtered by the processing and that of the input simulation. To simulate the cluster tSZ signal we used the universal pressure profile described in Arnaud et al. (2010) with the integrated tSZ signal 2 and redshift of CL J1226.9+3332 given in the Planck Collaboration XVIII. (2015) catalog, Y Planck 500 ∼ 3.82 × 10 −4 arcmin 2 and z = 0.89. We also add to the simulated cluster a Gaussian signal with flat spectrum (i.e. a random white noise) to explore angular scales at which the cluster signal is negligible. Up to now the NIKA2 SZ Large Program and NIKA analyses (Adam et al. 2015;Ruppin et al. 2018;Kéruzoré et al. 2020) considered one-dimensional transfer functions, hereafter 1D TF. In these cases circular symmetry is assumed and the 1D TF is obtained by averaging the power spectra ratio in Fourier-domain annuli at a fixed angular scale. Nevertheless, it is expected that the filtering is not isotropic in the map as it might depend on the scanning direction. This motivates the use of the two-dimensional transfer function, hereafter 2D TF. In the right panel of Fig. 4 we present the 2D TF describing the filtering in the NIKA2 150 GHz map of Fig. 2. The black line in the left panel of the figure shows the 1D TF, whereas the color lines correspond to the one-dimensional cuts of the 2D TF for the different directions represented in the right plot. Except for the scanning directions, the 2D TF is compatible with the 1D one and greater than 0.8 at large angular scales, meaning that the signal is well preserved. On the contrary, filtering is significant for 2 Y 500 is the integrated tSZ signal within a sphere of radius R 500 with D A the angular diameter distance at the cluster redshift. It is defined as:  angular frequencies below ∼ 0.5 arcmin −1 . At 0.4 arcmin −1 k 0.8 arcmin −1 the transfer function is larger than unity, meaning that the signal has been slightly enhanced by the data analysis process at these scales. In order to evaluate the impact of considering the anisotropy of the filtering, the analyses presented in the following sections will be carried out with both the 1D and 2D transfer functions (see Sect. 4.1.1).

Point source contamination
Point sources in the 150 GHz map contaminate the tSZ signal and also need to be considered. We start by identifying submillimetric sources by blindly searching for point sources in the NIKA2 260 GHz map of CL J1226.9+3332. By cross-checking the detections with a S/N greater than 3 with Herschel SPIRE 3 and PACS 4 catalogs, seven submillimetric sources were identified in the region covered by the NIKA2 maps. The position and fluxes from the above-mentioned catalogs for each submillimetric point source (PS1 to PS5, PS7 and PS8) are summarized in Table 1 and the corresponding Herschel names are given in Table 2.
For each of these sources, a modified blackbody spectrum model is adjusted to the fluxes in Table 1 together with the measurement of the flux at 260 GHz from the NIKA2 map. The spectra are then extrapolated to 150 GHz to obtain an estimate of the flux of each source at 150 GHz. In Table 2 we summarize the fluxes at 260 GHz obtained from the NIKA2 map and the extrapolated values at 150 GHz. The contribution of the noise and the filtering in the 260 GHz map are considered when measuring the fluxes. The estimates obtained with the AO and TO noise approaches give compatible results for all point sources. A more detailed explanation of the method used to deal with submillimetric point sources is given in Kéruzoré et al. (2020). The obtained probability distributions of the fluxes at 150 GHz are used as priors for the joint fit of the cluster pressure profile and the point sources fluxes described in Sect. 4.1.1.
Comparing results in Table 2, we notice very large uncertainties for the extrapolated flux of PS1 at 150 GHz when estimated with the TO noise map, but this will not affect substantially the following results for the simultaneous fit of the PS1 flux (see Table B.1) when estimating the cluster pressure profile. Moreover, we can compare its flux at 260 GHz to the measurement in Adam et al. (2015), where the source is called PS260. Using NIKA data, they obtain F 260 GHz = 6.8 ± 0.7 (stat.) ±1.0 (cal.) mJy which is consistent with our estimates.
Regarding PS6, it doesn't have a counterpart in Herschel SPIRE and PACS catalogs. But it appears as a weak signal in the Herschel maps, as a 3σ detection in the 260 GHz NIKA2 map and, moreover, it compensates the extended tSZ signal at 150 GHz (also clearly observed in Adam et al. (2015)). For this source the modified blackbody is used to get a prior knowlegde of the flux at 150 GHz for the assumed prior distributions of the spectral index and temperature (Kéruzoré et al. 2020). Another tricky point source is PS7. The extrapolated 150 GHz values shown in Table 2 are clearly overestimating the flux of the source. This is understandable since we do not have enough constraints for the low frequency slope of the spectral energy distribution. We choose to use the obtained values as upper limits of a flat prior for the flux of PS7 in the estimation of the cluster pressure profile.
In addition to submillimetric sources, according to the VLA FIRST Survey catalog (White et al. 1997), a radio source of 3.60 ± 0.13 mJy at 1.4 GHz is present in (R.A., Dec.) J2000 = (12h26m58.19s, +33d32m48.61s), hereafter PS9. This galaxy corresponds to the BCG identified in Holden et al. (2009) and the compact radio source detected with LOFAR in Di Gennaro et al. (2020). We know beforehand that the contribution of this radio source at 150 GHz is small, but given its central position it is important to consider it. Assuming a synchrotron spectrum F(ν) = F 0 (ν/ν 0 ) α with α = −0.7 ± 0.2, which describes the spectral energy distribution for an average radio source (Condon 1984), we get at 150 GHz, 0.1 ± 0.2 mJy. The obtained probability distribution of the flux at 150 GHz is also used as a prior for the fit in Sect. 4.1.1. The extrapolation of fluxes from radio to millimeter wavelengths can be dangerous and lead to biasing the electron pressure reconstruction. However, this is not the case for our analyses (see results in Sect. 4.1.1 and Table B.2). Notes.

XMM-Newton
Regarding X-ray data, CL J1226.9+3332 was observed by XMM-Newton (Obs ID 0200340101) for a total observation time of 90/74 ks (MOS/pn), reducing to 63/47 ks after cleaning. Data were reduced following standard procedures described in Bartalucci et al. (2017). The raw data were processed using the XMM-Newton standard pipeline Science Analysis System (SAS version 16). Only standard events from the EMOS1, 2, and EPN detectors with PATTERN < 4 and < 13, respectively, were kept. The data were further filtered for badtime events and solar flares. Vignetting was accounted for following the weighting scheme as described in Arnaud et al. (2001). Point sources were detected on the basis of wavelet filtered images in low ([0.3-2] keV) and high ([2-5] keV) energy bands, then subsequently masked from the events list. This process was controlled by a visual check allowing also to further extract obvious (sub-)structures present in the field. The instrumental background was modelled through the use of stacked filter-wheel closed observations, whilst the astrophysical contamination due to the Galaxy and the cosmic X-ray background were accounted for as a constant background in the 1D surface brightness analysis. This component was modelled in the spectral analysis following the method outlined in Pratt et al. (2010), using an annulus external to the target of radii 300"< θ < 480".

Electron pressure reconstruction from tSZ
The spectral distortion of the CMB caused by the thermal energy in the cluster, i.e. the tSZ effect, is characterized by its amplitude or Compton parameter, y (Sunyaev & Zeldovich 1972). This is directly proportional to the thermal pressure of the electrons in the ICM, P e , integrated along the line of sight, where σ T , m e and c are the Thomson cross section, the electron rest mass and the speed of light, respectively. Hence, the tSZ surface brightness is proportional to the Compton parameter integrated over the tSZ spectrum convolved by the NIKA2 bandpass and therefore, proportional to the integrated thermal pressure of the ICM in the cluster.

Pressure profile reconstruction with NIKA2
Reconstruction procedure To reconstruct the electron pressure in the ICM of CL J1226.9+3332 we fit a model map of the surface brightness of the cluster to the NIKA2 150 GHz map. The model map is obtained from the pressure profile of the galaxy cluster integrated along the line of sight in Compton parameter (y) units, following Eq. 1. We describe the pressure of the galaxy cluster with a Radially-binned spherical model Article number, page 7 of 23 A&A proofs: manuscript no. aanda (also known as Non-parametric model in Ruppin et al. (2017); Romero et al. (2018)): where P i and α i are the values of the pressure and the slope at the radial bin r i . The slope is directly calculated as: We initialize the pressure bin values taking random values from a normal distribution centered at the corresponding pressure from the universal profile of Arnaud et al. (2010) at each radial bin. The radial bins are chosen to cover mainly the range between the NIKA2 resolution and field of view capabilities. We center the pressure profile at the coordinates of the X-ray peak as determined from XMM-Newton data analysis (Sect. 3.1.2 and 3.2). The derived y-map is convolved with the NIKA2 beam, which is approximated by a two-dimensional Gaussian with FWHM = 17.6 (Perotto et al. 2020). In order to account for the attenuation or filtering effects due to data processing in the NIKA2 150 GHz map, the model map is also convolved with the transfer function (Sect. 3.1.4). We repeat this procedure for the 1D and 2D transfer functions. Finally, the y-map is converted into surface brightness units with a conversion coefficient, accounting for the tSZ spectra shape convolved by the NIKA2 bandpass, that is also left as a parameter of the fit (as done in Kéruzoré et al. 2020). Furthermore, for the comparison with the 150 GHz NIKA2 map, we added the contribution of point sources to the model map. Point sources are modeled as two-dimensional Gaussian functions and we repeat the procedure detailed in Kéruzoré et al. (2020) to fit the flux of each source at 150 GHz. Priors on the flux of the sources at 150 GHz are obtained from the results of the spectral fits presented in Sect. 3.1.5. The last component in the model map is a constant zero-level that we also adjust as a nuisance parameter.
The parameters (ϑ) of our fit are: the pressure bins describing the ICM of the cluster, the fluxes of the contaminant point sources, the conversion factor from Compton to surface brightness units and the zero-level. The likelihood that we use to compare our model M pixel by pixel to the data D is given by: Here C pix−pix is the pixel-pixel noise covariance matrix accounting for the residual noise in the NIKA2 150 GHz map (Sect. 3.1.3). We repeat the fit with the covariance matrices from both "angle order" and "time order" noise estimates.
We also compute the Y 500 integrated Compton parameter and compare it in the likelihood to the integrated Compton parameter measured by Planck Collaboration XVIII. (2015), Y Planck 500 = (3.82±0.79)×10 −4 arcmin 2 within an aperture of θ 500 = 1.907 arcmin. We do not compare the integrated Compton parameter at 5θ 500 as measured by Planck because it would require extrapolating the pressure profile far beyond the NIKA2 data.
For the map fit we use the PANCO2 pipeline (Kéruzoré et al. 2022) and follow the procedure described in Adam et al. (2015); Ruppin et al. (2018) and Kéruzoré et al. (2020). This pipeline performs a Markov Chain Monte Carlo (MCMC) fit using the emcee python package (Foreman-Mackey et al. 2019;Goodman & Weare 2010). The sampling is performed using 40 walkers and 10 5 steps, with a burn-in of 10 3 samples, and convergence is monitored following theR test of Gelman & Rubin (1992) and chains autocorrelation. The PANCO2 code has been successfully tested on simulations.

NIKA2 pressure radial profile
In order to estimate the robustness of the results of the above procedure, we have performed the fit to the NIKA2 data in four different cases with respect to the choice of noise residuals and transfer function estimates. Thus, we consider, AO1D (TO1D) and AO2D (TO2D) using the angle (time) order noise residual map and the 1D and 2D transfer functions, respectively. In Fig. 5 we compare the NIKA2 150 GHz map of CL J1226.9+3332 to the obtained best fit models and their residuals for these four analyses. Comparing the power spectra of the residual maps to the power spectra of the noise estimate maps, we see in Fig. 3 that for the TO case the fit residuals and the noise estimates power spectra are consistent. However, for the AO cases there is an excess of power in the fit residuals, which could be interpreted as coming from the signal due to the differential filtering effects that are not captured in the AO noise as discussed in Sect. 3.1.3. Regarding point sources, the reconstructed fluxes are consistent for the four analyses (see Tables B.1

and B.2).
We present in Fig. 6 the Radially-binned best-fit pressure profiles obtained for the four tested cases. The blue and cyan (dark and light green) dots correspond to the AO (TO) 1D and 2D transfer function estimates, respectively. The plotted uncertainties correspond to 1sigma of the posterior distributions derived from the MCMC chains. Overall, the four NIKA2 analyses give consistent results, specially in the radial ranges in which we expect the NIKA2 results to be reliable, i.e. between the beam and the field of view (FoV) scales, both represented with dashed vertical lines in the figure. We give the HWHM of the NIKA2 beam (17.6"/2) and half the diameter of the FoV (6.5'/2) in the physical distances corresponding to the redshift of the cluster.
In terms of noise estimates we observe that the uncertainties on the pressure bin estimates are slightly larger for the time ordered case as one would expect. However, we notice no significant bias between the time and angle ordered results. The effect of the transfer function is hard to evaluate: even if the 2D TF is a more precise description of the filtering in the map, when fitting a spherical cluster model the use of the 1D TF gives consistent results. In the following we use the results for the four analyses to evaluate possible systematic uncertainties induced by the NIKA2 processing. Fig. 6, we compare our results to the profiles obtained in R18 with tSZ data from NIKA (pink), Bolocam (black) and MUSTANG (yellow). MUSTANG's high angular resolution (9" FWHM at 90 GHz) enables to map the core of the cluster, whereas Bolocam's large field of view (8' at 140 GHz) allows one to recover the large angular scales. NIKA and the improved NIKA2 camera are able to cover all the intermediate radii. The consistency of the different pressure bins in the radial range from the NIKA2 beam to the FoV proves the reliability of the reconstruction with NIKA2 data.  Before going further, we have to consider again the effect of the filtering on the NIKA2 data. The filtering due to the data processing affects mainly small angular frequencies, i.e. small k numbers (Sect. 3.1.4), which is translated into large angular scales in real space. In this case it means that the region at ∼ 1000 kpc from the center of the cluster is strongly filtered. For this reason, we cast doubt on the results of our fits for the last NIKA2 bin in pressure.

Thermodynamical profiles from X-rays
The electron density and temperature profiles were extracted following the methodology described by Pratt et al. (2010) and Bartalucci et al. (2017). In short, the vignetted-corrected and background-subtracted surface brightness profile obtained in concentric annuli from the X-ray peak is deconvolved from the PSF and geometrically deprojected assuming spherical symmetry using the regularisation technique described in Croston et al. (2006).
The temperature profile is derived in bins defined from the afore-derived binned surface brightness profile through a Article number, page 9 of 23 A&A proofs: manuscript no. aanda spectral analysis modelling the ICM emission via an absorbed MEKAL model under XSPEC 5 and accounting for both the instrumental and astrophysical backgrounds. The derived 2D temperature profile is then PSF-corrected and deprojected following the "Non-parametric-like" method presented in Bartalucci et al. (2018).
The gas pressure, P, and entropy, K, profiles were then derived from the deprojected density, n e , and temperature, T , profiles assuming P = n e × T and K = T/n 2/3 e , respectively. In this paper we focus on the former, which is shown in Fig. 6. as open black circles. The electron density and temperature profiles are shown in Fig. 7.

Hydrostatic equilibrium
Under the hydrostatic equilibrium (HSE) hypothesis we can compute, for a spherical cluster, the total cluster mass enclosed within the radius r as: M HSE (< r) = − 1 µm p G r 2 n e (r) dP e (r) dr (5) where µ, m p and G are the mean molecular weight of the ICM gas, the proton mass and the gravitational constant, respectively. We assume µ ≈ 0.6 (Pessah & Chakraborty 2013; Ettori et al. 2019) for the gas. Combining the pressure profiles obtained from the thermal SZ or X-ray data with the electron density from the X-ray we can reconstruct the mass of the galaxy cluster as in

Pressure profile modeling for mass estimation
Deriving the mass directly from the Radially-binned profiles presented above (see Fig. 6) leads to non-physical results (i.e. negative mass contributions) as no condition was imposed regarding the slope of the profile in the pressure reconstruction. This was done to prevent extra constraints on the pressure profile induced by assumptions on the model. To overcome that issue, we fit here pressure models ensuring physical mass profiles to the Radiallybinned tSZ results from Sect. 4.1.1. We will consider two differ-5 https://heasarc.gsfc.nasa.gov/xanadu/xspec/ ent approaches: 1) a direct fit of a generalized Navarro-Frenk-White (gNFW) profile to the Radially-binned pressure, and 2) an indirect fit of a Navarro-Frenk-White (NFW) mass density model under the HSE assumption. In both cases and aiming for a precise reconstruction of the HSE mass, which requires having accurately constrained slopes for the pressure profile, we combine the NIKA2 pressure bins with the results obtained in R18 6 .

gNFW pressure model
The first approach, which has been used in previous NIKA2 studies Kéruzoré et al. 2020;Ferragamo et al. 2022), consists in fitting the widely used gNFW pressure profile model (Nagai et al. 2007) to the tSZ data: where P 0 is the normalization constant, b and c the external and internal slopes, r p the characteristic radius of slope change and a the parameter describing the steepness of the slopes transition. We perform a MCMC fit using emcee package. The likelihood function is given by: where P N2 and C represent the NIKA2 Radially-binned pressure profile bins and associated covariance matrix. P R18 k and ∆P R18 k are the R18 pressure profile data and uncertainties. And M. Muñoz-Echeverría et al.: Multi-probe analysis of the galaxy cluster CL J1226.9+3332 Fig. 8. Pressure profile and best-fit for the gNFW (left) and NFW (right) models. The data points correspond to the NIKA2 Radially-binned results for the four data sets discussed above, and to the NIKA, MUSTANG and Bolocam bins from R18. Blue and green solid lines represent the best-fit values for the four NIKA2 pressure estimates considered. The shaded regions show the 2.5th, 16th, 84th and 97.5th percentiles. P gNFW (ϑ) are the gNFW pressure profile values for a set of parameters ϑ = [P 0 , r p , a, b, c]. As we do not rely on the value of the last NIKA2 pressure bin, we have chosen to modify the NIKA2 inverse covariance matrix C −1 by setting the last diagonal term to [C −1 ] 6,6 = 0, so that the correlation of the last bin with the others is taken into account, but not its value. We also set a constraint on the integrated Compton parameter of the model Y gNFW 500 (ϑ) using again the Planck satellite PSZ2 catalog results, Y Planck 500 (Planck Collaboration XVIII. 2015). However, we find that the impact of this constraint is completely negligible for this cluster. Furthermore, we have added an extra condition to the fit to ensure increasing HSE mass profiles with radius: r 2 n e (r) dP e (r) dr < 0. The best-fit gNFW pressure profiles (solid lines) and uncertainties (shaded area) are presented in the left panel of Fig. 8 for the four sets of NIKA2 data discussed in Section 4.1.1. We observe that the best-fit models are a good representation of the data over the full range in radius as demonstrated by the corresponding reduced χ 2 , which are close to 1 for all the cases (see solid lines in Fig. C.2). The posterior distributions of the ϑ gNFW parameters can be found in Appendix C.
In addition, it is interesting to compare our results to those from Planck for which a similar modeling was used. In Fig. 9 we present the 2D posterior distributions of the integrated Compton parameter at 5R 500 (with R 500 calculated independently in each case) with respect to the Θ s parameter of the gNFW model, at 68%, 95% and 99% C.L.. Θ s and r p are related via the angular diameter distance at the cluster redshift: tan(Θ s ) = r p /D A . We compare the results obtained in Planck Collaboration XVIII. (2015) (with the MMF3 matched multi-filter, available in the Planck Legacy Archive 7 ) to the constraints from the gNFW profiles obtained in this work with NIKA2, R18 and XMM-Newton data. Our contours were obtained varying all the parameters in the gNFW model fit, while for Planck a, b and c were fixed. For simplicity, we only show the contours for the NIKA2 AO1D case. This figure illustrates the important gain in precision due to 7 https://pla.esac.esa.int/#catalogues Fig. 9. Distribution of Y 5R500 with respect to Θ s for the gNFW pressure model fits to Planck data (grey) in Planck Collaboration XVIII. (2015) and to the NIKA2 + R18 + XMM-Newton data (blue) in this work. Different contours show the 68%, 95% and 99% confidence intervals. The black star corresponds to the intersection between the Planck Collaboration XVIII. (2015) distribution and the X-ray scaling law shown in Fig. 16 in Planck Collaboration XVIII. (2015). high resolution observations: resolving the galaxy cluster allows us to determine, even at such high redshift, the Θ s characteristic radius. Contours are marginally in agreement.

NFW density model
We present here a different approach for the modeling of the pressure profile in the scope of estimating the HSE mass. Indeed, the estimation of the pressure derivative (Eq. 5) can be very problematic as: 1) it is very sensitive to local variations in the slope of the pressure profile and, 2) it requires, as discussed above, additional constraints to ensure recovering physical masses. To over-come these issues we model the pressure profile starting from a mass density model and assuming HSE. An equivalent idea is the "backward process" to fit X-ray temperatures described in Ettori et al. (2013) and references therein. This method was used for the mass reconstruction in Ettori et al. (2019) and in Eckert et al. (2022). From the HSE defined in Eq. 5, we can write: M HSE (< r) r 2 dr (8) Moreover, we can relate a radial mass density profile ρ(R) to the mass by, that allows us to relate the pressure directly to a mass density profile. We use here the NFW model, which is a good description of dark matter halos (Navarro et al. 1996) and has been widely used in the literature (e.g. Ettori et al. 2019): where ρ c is the critical density of the Universe at the cluster redshift and δ c 200 8 is a function that depends only on c 200 , the concentration parameter (we switch here from an overdensity of 500 to 200 in order to conform to most of previous works). Finally, r s represents a characteristic radius and it is also a free parameter of the model. Using this definition we obtain, P zero − P(r a ) = −µm p G4πρ c δ c 200 (c 200 )r 3 s r zero r a n e (r) r 2 1 1 + r/r s + ln(1 + r/r s ) − 1 dr (11) where r zero is a radius at which we are dominated by a zero level component. We perform a MCMC analysis similar to the one described above for the gNFW pressure profile model. In this case the free parameters of the model are ϑ = [c 200 , r s , P zero ]. At each step of the MCMC we compute the integral in Eq. 11 to evaluate P(ϑ) as needed for the likelihood function in Eq. 7. Calculating the integral can be computationally very expensive. As the result of this integral depends only on r s and r a , we create a grid of the integrals for a range of r s values (from 100 to 2000 kpc) and r a the radial bins of interest. We use this grid to interpolate the values of the integrals at each step. Flat priors are given for the concentration, 0 < c 200 < 4. We also make use of the python NFW package 9 .
The best-fit pressure profiles and uncertainties are presented in the right panel of Fig. 8 for the four NIKA2 Radially-binned data sets discussed above. The posterior probability distributions of the free parameters of the models are shown in Appendix C. The posterior distributions of the c 200 and r s parameters can be compared to the results for the analysis of clusters in X-rays in Pointecouteau et al. (2005); Ettori et al. (2019); Eckert et al. (2022). In these studies, c 200 spans from 1 to 6 and r s from 200 kpc to 1200 kpc, which is compatible with our results. We note that our constraints on c 200 do not significantly differ from our prior on this parameter. We find that the NFW model is overall a good fit to the data as shown by the reduced χ 2 ∼ 1 (see   10. HSE mass profile estimates for CL J1226.9+3332 obtained with the NIKA2 (angle order 1D) and R18 tSZ data combined with the XMM-Newton electron density profile. The solid magenta and dashed blue lines correspond to gNFW and NFW methods, respectively. The shaded areas show the 2.5th, 16th, 84th and 97.5th percentiles. Empty dots correspond to the HSE mass profile obtained from XMM-Newton data only. Vertical dashed lines show the R HSE 500 obtained from each mass profile. The grey region represents the radial ranges at which the profiles are extrapolated.
increase significantly in the outskirt of the cluster with respect to the gNFW pressure profile model discussed in the previous section. This can be probably explained by the flexibility of the NFW-based approach, which is high enough to show that the last point in the profile is not well-constrained by the data.

HSE mass estimates
We present in Fig. 10 the HSE mass profiles inferred from the gNFW best-fit pressure profile in combination with the XMM-Newton electron density, and, from the NFW density best-fit model. Uncertainties (shaded areas) are obtained directly from the MCMC chains by computing the HSE mass profile for each sample from the model parameters. For the sake of clarity we only present the masses obtained with NIKA2 AO1D estimates, but we change the color code for gNFW profile so that we can differentiate both results.
The capability of the pressure model to describe the shape of the profile slopes is the key element for a good mass reconstruction. From the above results it seems that the NFW approach does not have enough degrees of freedom to fully describe slope variations in the reconstructed pressure profile. Nevertheless, the resulting HSE mass profiles for both models are compatible within 2σ. The vertical dashed lines in the figure represent the R HSE 500 for each mass profile. Slight differences in the shape of the pressure profile at these radial ranges are critical for defining R HSE 500 . These mass profiles, obtained from the combination of tSZ and X-ray data, are also compared to the X-ray-only HSE mass estimate in Fig 10. Assuming spherical symmetry, the X-ray mass profile was derived, following the Monte Carlo procedure detailed in Démoclès et al. (2010) and Bartalucci et al. (2017) with the XMM-Newton electron density and temperature profiles presented in Sect. 4.2. Despite the different behaviour of the X-ray-only profile in the cluster core, it is consistent with the tSZ+X estimates at around R HSE 500 .
Article number, page 12 of 23 M. Muñoz-Echeverría et al.: Multi-probe analysis of the galaxy cluster CL J1226.9+3332 Fig. 11. 1D and 2D probability distributions for M HSE 500 and R HSE 500 from the gNFW (left) and NFW (center) models in the combined XMM-Newton and NIKA2 and R18 data, and, from the XMM-Newton X-ray data only (right). The different blue and green lines correspond to results for the four NIKA2 test cases considered. From the reconstructed HSE mass profiles we can obtain R HSE 500 − M HSE 500 probability distributions for each of the considered cases. We present in the left and central panel of Fig. 11 the R HSE 500 − M HSE 500 distributions for the gNFW and NFW models. They have been obtained from the MCMC chains in the same way uncertainties in Fig. 10 have been computed. The width of the ellipses is an artifact from the display procedure and each value of M HSE 500 is associated to a single value of R HSE 500 .
We present here results for the four NIKA2 analyses (AO1D/2D and TO1D/2D). The results are consistent, with small dependency on the chosen TF estimate. From the comparison of left and central panels in Fig. 11 we verify that the largest uncertainty in the HSE mass estimates comes from the modeling of the pressure profile. In spite of this effect, the reconstructed HSE mass profiles are compatible within 1σ. The right panel of Fig. 11 shows the R HSE 500 − M HSE 500 probability distribution obtained with XMM-Newton data only. Even if it is compatible with the gNFW and NFW results, the X-ray-only results prefer lower HSE masses. A similar effect was observed for ACT-CL J0215.4+0030 cluster (Kéruzoré et al. 2020), but not for PSZ2 G144.83+25.11 . We summarize in Table 3 the marginalized M HSE 500 masses obtained in this work. We give the mean value and the 84th and 16th percentiles. For gNFW and NFW we combine the probability distributions obtained for the four NIKA2 results so that the results account for the systematic effects from NIKA2 data processing.
In general, we conclude that minor variations in the HSE mass profiles at ∼ R 500 can lead to large variations on M 500 . We also stress that the degeneracy between M 500 and R 500 is intrinsic and not induced by systematic or stochastic uncertainties.

Lensing mass
Lensing masses, by contrast to HSE ones, probe the total mass without the assumptions on the dynamical state of the cluster. For this reason, it is of great interest to compare the HSE masses to lensing estimates. In this section we present the lensing mass estimate for CL J1226.9+3332.

Lensing data
We use the CLASH convergence maps (hereafter κ-maps) obtained from the weak and strong lensing analysis by Zitrin et al. (2015). In the analysis the authors reconstructed the κ-maps for the 25 massive CLASH clusters (Postman et al. 2012) using two different lensing models: 1) Light Traces Mass, LTM and 2) Pseudo Isothermal Elliptical Mass Distribution plus an elliptical NFW dark matter halo, PIEMD+eNFW. These models have been detailed in previous studies (Zitrin et al. 2009;Zitrin et al. 2013aZitrin et al. ,b, 2015. In this work we compute the mass estimate for both κ-maps in order to account for differences in the convergence map modelling.

Lensing mass density profile
For the lensing mass profile reconstruction we followed the approach described in Ferragamo et al. (2022). The convergence maps describe the projected mass density of the cluster, Σ, in critical density units, κ = Σ/Σ crit , with Here D s , D l and D ls correspond to the angular diameter distance between the observer and the source, the observer and the lens (the cluster) and between the source and the lens, respectively. The publicly available CLASH κ-maps (Zitrin et al. 2015) have been normalised to D s /D ls = 1.
To estimate the lensing mass profile of CL J1226.9+3332, we fit a mass density model to the Σ-map. We assume spherical symmetry and a NFW density profile. We choose to directly fit the analytical projected NFW density profile (Eq. 5 in Ferragamo et al. (2022) derived from Navarro et al. (1996); Bartelmann (1996)) to the radially averaged projected profiles of the Σ-maps. We consider as free parameters r s and c 200 (see Eq. 10). The fit is performed via a MCMC analysis using the emcee software Fig. 12. Projected mass density profiles obtained from CLASH convergence maps for the LTM (left) and PIEMD+eNFW (right) models. We also show the best-fit NFW model (orange and red lines) and the 2.5th, 16th, 84th and 97.5th percentiles (shaded area). and the NFW python package. We center the projected mass profiles at the same position as for the pressure and electron density profiles in Sect. 4.1 and 4.2, the X-ray center. Uncertainties were computed from the dispersion in each radial bin and for this analysis we have not considered the possible correlation between the bins.
We show in Fig. 12 the radial profiles of the projected mass density for CL 1226.9+3332 as obtained from the CLASH LTM (left) and PIEMD+eNFW (right) convergence maps. We present for both profiles the projected best-fit NFW density model and percentiles (shaded area). We observe that for the LTM convergence map the best-fit NFW model underestimates the data except for cluster core and that the uncertainties in the model do not fully account for this. By contrast, the fit for the PIEMD+eNFW succeeds in representing the data.
From these results we conclude that the main uncertainties in the reconstruction of the lensing density profile comes from the reconstruction of the convergence map. Thus, in the following we will account for those.

Lensing mass estimates
From the obtained NFW density profiles we can reconstruct the lensing mass profiles (Eq. 9) and subsequently the M lens 500 and R lens 500 probability distributions. We present in Fig. 13 the main results of this analysis for both convergence map models: LTM and PIEMD+eNFW. We obtain consistent results and accounting for their combined posterior distributions we measure M lens 500 = 7.35 +0.53 −0.51 × 10 14 M mean value and the 84th and 16th percentiles.
The Jee & Tyson (2009) weak-lensing analysis did not provide the direct M lens 500 , but evaluated the lensing mass at the R 500 from Maughan et al. (2007). This result, shown in purple in Fig. 13, is consistent within uncertainties with our results when evaluating the lensing mass at the same radius for the LTM and PIEMD+eNFW analyses (see orange and red error bars in the figure). Merten et al. (2015) performed an independent analysis of the CLASH data, reconstructing their own convergence map. The projected mass density profile presented in Fig. 16 in Merten et al. (2015) shows a denser cluster than the profiles from the convergence maps used in this work (Fig. 12). For this rea- Fig. 13. Probabiliy distribution for M lens 500 and R lens 500 obtained from the fit of the NFW density profile model on the CLASH PIEMD+eNFW (red) and LTM (orange) convergence maps. We show 1 and 2σ C.L. contours. The purple cross corresponds to the lensing mass estimate in Jee & Tyson (2009). We also show for comparison lensing masses estimated at 880 kpc for the PIEMD+eNFW (red) and LTM (orange) as vertical error bars at 1σ C.L. The brown stars correspond to M lens 500 from literature. We differentiate the mass from Merten et al. (2015) (full) and Sereno (2015); Sereno & Covone (2013) (empty). son, Merten et al. (2015) obtained, also with a NFW density fit, 35% larger masses than in Jee & Tyson (2009). The corresponding M lens 500 is shown with a full brown star in Fig. 13. The disturbed state of CL J1226.9+3332 could be the reason, according to Merten et al. (2015), for the different lensing mass estimates. Moreover, the high redshift of the cluster makes more difficult the precise reconstruction of the convergence map. The empty brown star in Fig. 13 corresponds to the lensing mass from Co-MaLit (Sereno 2015) obtained from the analysis in Sereno & Covone (2013). In this case, the result is compatible with our mass estimates.

Comparison of mass estimates
7.1. CL J1226.9+3332 mass at R 500 The comparison of different mass estimates is difficult and can lead to wrong physical conclusions. In particular, when comparing integrated masses the radius at which the mass is computed has a significant impact: R 500 and M 500 being constrained at the same time, we are affected by a degeneracy. For this reason, in Fig. 14 we show, in the R 500 − M 500 plane, the results from the literature with the ones obtained in this work.
The green and blue contours show the R HSE 500 − M HSE 500 results obtained in this work for the gNFW (solid lines) and NFW (dashed lines) tSZ and X-ray data combined analyses. For comparison the full grey diamonds correspond to the results from the literature presented in Fig. 1 also for combined tSZ and X-ray data. We observe that the results in this paper are compatible with previous analyses within 1σ, centered around ∼ 7 × 10 14 M . Regarding X-ray-only results, the HSE mass estimates obtained in this work with XMM-Newton data (grey contours) suggest mass values centered at ∼ 5 × 10 14 M . This is in agreement with the lowest estimates from the literature (open grey diamonds) presented in Bulbul et al. (2010) and Maughan et al. (2007). On the contrary, the results from Mantz et al. (2010) and Mroczkowski et al. (2009) show larger masses. However, the M HSE 500 in Mantz et al. (2010) is not a direct measurement, but an extrapolation from a gas mass measured at R 2500 converted into total mass, making this result less reliable. Overall, for CL J1226.9+3332 the HSE masses obtained only from X-ray data tend to lower values than those from the combination of tSZ and X-rays.
The result from Planck Collaboration XVIII. (2015) (lowest full diamond) is also a special case, as it is not a direct mass measurement, but a mass obtained from the X-ray-derived scaling relation (Eq. 7 in Planck Collaboration XX. (2013)) applied to the SZ measurement. This may explain why it lies at the border between the X-ray-only data and the tSZ+X combined results. The differences observed between X-ray-only and the combined tSZ+X results could have a physical and observational origin. For such a high redshift cluster X-ray observations become challenging. If the southwestern sub-clump in the cluster is really a hot but not dense structure (as suggested by Jee & Tyson (2009)), the electron density measurements from X-ray observations might be a struggle.
We also show in Fig. 14 the lensing, dynamical and Virial mass estimates. The lensing mass estimates from this work for the CLASH LTM and PIEMD+eNFW convergence maps are presented as dark orange and red contours, respectively. We observe that they are consistent with the lensing mass from Sereno (2015); Sereno & Covone (2013) as well as with HSE mass estimates, but very different from the Merten et al. (2015) lensing estimate for the reasons explained in Sect. 6. The Virial masses estimated in Mroczkowski (2011) and Mroczkowski (2012) without X-ray data are shown as pink squares. They rely on the Virial relation and on given pressure and density profile models to relate directly the integrated tSZ flux to the mass (Eq. 15 in Mroczkowski (2011)). This kind of analysis seems a good alternative to the HSE mass for clusters without X-ray data. The dynamical mass estimates (purple circles), which we would expect to be larger than the HSE estimate, appear particularly low for CL J1226.9+3332 (Aguado- Barahona et al. 2022 (13) where M true is the total true mass of the cluster. From the observational point of view there are hints of a non null hydrostatic bias. As mentioned in Sect. 1, one example is the tension observed between the cosmological parameters derived from Planck cluster number counts and those from the CMB analyses (Planck Collaboration XXI. 2013). A possible explanation for this tension is that cluster masses are underestimated. According to Planck Collaboration XXIV. (2015) the bias needed to reconcile the cosmological constraints obtained from the CMB power spectra to the cluster counts is 1 − b = M HSE 500 /M true 500 = 0.58 ± 0.04. A compatible value was obtained from the updated analysis in Salvati et al. (2019), 1 − b = M HSE 500 /M true 500 = 0.62 ± 0.05. The hydrostatic bias has also been the topic of a large number of studies on numerical simulations (see Ansarifard et al. 2020;Gianfagna et al. 2021, and references therein). However, simulation-based analyses agree on values in the range of 0.75 to 0.9 for 1 − b, not able, therefore, to reconcile the mentioned tension.
In a hierarchical formation scenario one would expect clusters at higher redshift to be more disturbed and therefore the hydrostatic equilibrium hypothesis to be less valid (Neto et al. 2007;Angelinelli et al. 2020). Thus, this possible redshift evolution has been studied in the literature: e.g. Salvati et al. (2019) find a modest hint of redshift dependence for the bias. However, other works based on cluster observations (McDonald et al. 2017) do not find traces of evolution of the morphological state with redshift. To date, the possible bias dependence with redshift is not confirmed in simulations (see Gianfagna et al. 2021, and references therein).

Hydrostatic-to-lensing mass bias estimates
From the observational side the real HSE bias is unachievable, as one cannot access the true mass of a cluster. But it can be approximated using mass estimates that do not rely on the HSE hypothesis and trace the total mass of the cluster, for instance the lensing mass. In this work we have computed the hydrostaticto-lensing mass bias using the results obtained in Sect. 5.3 and 6 (see Sereno & Ettori (2015) for an analysis of the CoMaLit samples). For the lensing mass M lens 500 we combine the probability distributions of both lensing models in Fig. 11. On the contrary, we consider the different HSE mass estimates independently. Assuming that HSE and lensing masses are uncorrelated estimates, we have combined their probability distributions and computed the ratio, M HSE 500 /M lens 500 = 1 − b HSE/lens . Fig. 14. R 500 − M 500 plane summarizing the results for CL J1226.9+3332. In the case of the HSE mass estimates the blue and green 1σ contours show the results obtained in this work combining tSZ and X-ray data for the four NIKA2 analyses. The solid and dashed lines are for the gNFW pressure and the NFW density models, respectively. The grey contour corresponds to the HSE mass estimate for the XMM-Newton data-only. In the case of literature data the full (empty) grey diamonds represent HSE masses from the combination of tSZ and X-ray data (X-ray-only results). The red and orange contours correspond to the lensing mass estimates obtained from the CLASH LTM and PIEMD+eNFW convergence maps in this work, respectively. Pink squares show the tSZ-only mass assuming Virial relation, purple circles are dynamical mass estimates and the brown stars the lensing estimates. The diagonal bright grey line defines the R 500 − M 500 relation. Slight deviations from this line are only due to differences in the cosmological model used in each work. We present in Fig. 15 the hydrostatic-to-lensing mass ratio at R 500 . The same color code as in previous figures is used to distinguish the HSE estimates that have been obtained with each of the NIKA2 noise and filtering estimators. We draw with solid lines the results obtained modeling the pressure profile with the gNFW model and with dashed lines the results from the NFW fit. All of them are inferred from the NIKA2, R18 and XMM-Newton data. The grey contours show the bias using the X-rayonly HSE mass results. We present in Table 4 the marginalized hydrostatic-to-lensing mass ratio and uncertainties for the different cases considered. The results for the gNFW and NFW cases correspond to a combination of the four NIKA2 analyses. As concluded in Ferragamo et al. (2022), we observe that the value of the bias is sensitive to the considered data sets and modeling choices.

Summary and conclusions
The precise estimation of the mass of single clusters appears extremely complicated and affected by multiple systematic effects (Pratt et al. 2019), but it is the key for building accurate scaling relations for cosmology. The NIKA2 SZ Large Program aims at providing a SZ-mass scaling relation from the combination of NIKA2 and XMM-Newton data. Within this program, we presented a thorough study on the mass of the CL J1226.9+3332 galaxy cluster.
We obtained NIKA2 150 and 260 GHz maps, that allowed us to reconstruct the Radially-binned pressure profile for the ICM from the tSZ data. To characterize the impact of the data processing, we repeated the whole analysis for two pipeline-filtering transfer functions and noise estimates for the 150 GHz map. We accounted for the presence of point sources that contaminate the negative tSZ signal at 150 GHz. The reconstructed NIKA2 pressure bins are compatible, within the angular scales accessible to NIKA2, with the profiles obtained from three independent instruments in R18. This validates the pressure reconstruction procedure that will be used for the whole sample analysis in the NIKA2 SZ Large Program.
We compared two approaches to estimate the HSE mass from the combination of tSZ-obtained pressure and X-ray electron density profiles. We considered either a gNFW pressure model (traditionally used for this kind of analyses) or an integrated NFW density model. The second seems a promising ap- Fig. 15. HSE mass estimates with respect to lensing ones at R 500 . Blue and green solid lines correspond to the HSE masses obtained from the gNFW pressure model fit and dashed lines to the NFW method. Grey line correspons to the HSE mass obtained with XMM-Newton only data. We present the 1σ and 2σ C.L. contours. The lensing mass distribution is the combination of the results for the PIEMD+eNFW and LTM analyses.
proach to ensure radially increasing HSE mass estimates. However, other density models should be tested, in order to describe more satisfactorily the shape of the pressure profile. Both methods give completely compatible HSE mass profiles and integrated M HSE 500 . From the comparison of the different mass estimates, we also conclude that for the moment, when estimating the HSE mass in the NIKA2 SZ Large Program, the error budget is dominated by model dependence rather than by the instrumental and data processing systematic effects that we investigated. In addition, these results are in agreement with the X-ray-only HSE mass estimate obtained in this paper from the XMM-Newton electron density and temperature profiles. Nevertheless, the latter prefers lower mass values than the combined tSZ+X-ray results.
We have also shown that our results are compatible with all the HSE mass estimates found in the literature within uncertainties, which are large. We think that the only way to reduce the current uncertainties is to precisely constrain the slope of the mass profile at ∼ R 500 , as we have proved that very similar mass profiles overall can result in significant differences at M 500 .
From two differently modeled CLASH convergence maps, we have reconstructed the lensing mass profile of CL J1226.9+3332 and measured the hydrostatic-to-lensing mass bias. We have found that for CL J1226.9+3332 the bias is b HSE/lens ∼ 0.1 for tSZ+X-ray combined HSE masses and ∼ 0.4 for X-ray-only estimates, with ∼ 0.1 uncertainties depending on the models. In spite of the large uncertainties, we have the sensitivity to measure the bias of a single cluster, even for the highest-redshift cluster of the NIKA2 SZ Large Program. This is the second cluster (the first one was studied in Ferragamo et al. (2022)) in the NIKA2 SZ Large Program with such a measurement and the analysis of the HSE-to-lensing bias for a larger sample will be the topic of a forthcoming work (as done in Bartalucci et al. (2018) with X-ray data). Measuring the hydrostatic-to-lensing mass bias associated with LPSZ clusters will bring a new perspective to the bias from the tSZ side for spatially resolved clusters at intermediate and high redshifts.  1500.0 2 ± 2 4 ± 2 4 ± 3 8 ± 3