© The Authors 2024

1. Introduction

The distribution of galaxy clusters in mass and redshift is a key tool for estimating cosmological parameters (Vikhlinin et al. 2009; Planck Collaboration XX 2014; Costanzi et al. 2019). Nevertheless, the mass of clusters is not an observable quantity and has to be estimated under several hypotheses from observations.

Firstly, some methods assume that galaxy clusters are in hydrostatic equilibrium (HSE) and combine intra-cluster medium (ICM) observables, such as the X-ray emission and/or the Sunyaev-Zel’dovich (SZ) effect (Sunyaev & Zeldovich 1972; Planck Collaboration XI 2011), to infer the HSE mass (Mroczkowski et al. 2009; Adam et al. 2015; Ruppin et al. 2018; Kéruzoré et al. 2020; Eckert et al. 2022a; Muñoz-Echeverría et al. 2023). Secondly, others relate the velocity dispersion of the galaxies in the cluster to its total mass (Aguado-Barahona et al. 2022; Biviano & Girardi 2003). Thirdly, the lensing effect on background sources due to the gravitational field of the cluster can be used to reconstruct the mass (Merten et al. 2015; Zitrin et al. 2013). Each method is known to bias the results differently (Pratt et al. 2019) due to instrumental limits or departures from the assumed hypotheses, and comparison works making use of different observational results try to understand those issues (e.g. Sereno & Ettori 2015).

However, other biases are intrinsic to cluster physics and geometry and need to be quantified from simulations. Simulations allow one to compare the three-dimensional properties of clusters, needed for cosmological studies, to those inferred from projected maps.

The orientation of clusters is known to affect the mass reconstruction. As shown in Oguri et al. (2009), Oguri & Blandford (2009), Hennawi et al. (2007), Meneghetti et al. (2010a), Gralla et al. (2011), and Euclid Collaboration (2024), when clusters are elongated along the line of sight, masses are overestimated. The contrary happens when the major axes of clusters are on the plane of the sky. Thus, the projections through which clusters are observed impact their mass estimates. In Meneghetti et al. (2014), the authors explored the impact of projection effects by using the MUSIC-2 simulations (Sembolini et al. 2013; Biffi et al. 2014), that are N-body MultiDark simulations¹ where baryons were added using smoothed-particle hydrodynamics (SPH) techniques on the dark matter, in the non-radiative flavour and without accounting for the energy feedback from active galactive nuclei (AGN). The authors estimated the mass by fitting different density models to the spherical mass density profiles of the clusters as well as to the projected mass maps. Overall, masses reconstructed from projected maps are ∼13 − 14% more scattered than those estimated from the fit of three-dimensional density profiles and, according to this work, masses estimated from projected data are on average under-estimated by 5%. They affirmed, as in Giocoli et al. (2012a), that this bias is due to the orientation of clusters, which would be preferably elongated on the plane of the sky. Lensing observations give an estimate of the projected matter density distribution; therefore, the difficulty resides in recovering precise three-dimensional profiles.

The prolateness was also presented in Giocoli et al. (2012a) as a source of bias. According to simulations in the cold dark matter (CDM) framework (Shaw et al. 2006), clusters are more frequently prolate systems than oblate-shaped. In addition, the presence of substructures was found in Giocoli et al. (2012a) as the second contributor to the mass bias, after halo triaxiality.

The model chosen for the mass reconstruction can also be a source of bias. In Meneghetti et al. (2014), the authors used the Navarro-Frenk-White (NFW, Navarro et al. 1996), generalised Navarro-Frenk-White (gNFW, Nagai et al. 2007; Zhao 1996), and Einasto (Einasto 1965) profile models to fit mass density radial profiles in ranges between 0.02R_vir and R₂₀₀. We note that R_vir is the virial radius of the cluster (Peebles 1980; Bryan & Norman 1998) and R₂₀₀ is the radius at which the mean density of the cluster is 200 times the critical density of the Universe².

They concluded that, as expected, models with three parameters (gNFW and Einasto) fit the density profile better than those with two (NFW). An important projection effect was also found regarding modelling: not all clusters that have density profiles following a NFW shape in three dimensions have a NFW-like density when projected (Meneghetti et al. 2014). In Euclid Collaboration (2024), the authors created Euclid-like weak lensing observables from THE THREE HUNDRED simulation shear maps (see Sect. 2.2.1 for more details on the simulations). By fitting smoothly truncated NFW (tNFW) density models (Baltz et al. 2009) to these data, they concluded that the chosen truncation radius in the model impacts the bias of reconstructed masses. In addition, they showed that the bias at R₂₀₀ and its relative uncertainty are smaller if the concentration parameter is fixed, namely at c₂₀₀ = 3.

Moreover, as presented in Rasia et al. (2012), another reason to explain the bias is the approach used to select the sample. For clusters selected according to their X-ray luminosity, the reconstructed concentration-mass relation has a larger normalisation and steeper slope than the genuine relation (Meneghetti et al. 2014; Rasia et al. 2013). The reason behind this is that, for a given mass, the most luminous clusters are the most concentrated ones. On the contrary, if clusters are chosen for their strong lensing signal, they are preferentially elongated along the line of sight (Meneghetti et al. 2010a; Giocoli et al. 2014) and, as a consequence, the masses obtained from projected maps are overestimated.

Many are the effects at the origin of the bias and the scatter of the weak lensing masses of clusters. As summarised in Lee et al. (2023) and references therein, assuming a NFW profile and a concentration-mass relation can introduce uncertainties in the mass from 10% to 50% due to non-sphericity of clusters, structures along the line of sight, miscentring and halo concentration. Regarding the concentration-mass relation, from an analysis of mock galaxy cluster lenses created with the MOKA³ code (Giocoli et al. 2012a,b) concluded that the amplitude of the relation is lower when derived from lensing analyses, therefore, from projected mass tracers, than from the three-dimensional studies.

In Rasia et al. (2012), 20 haloes simulated to mimic Subaru observations (Fabjan et al. 2011) were used with three projections per cluster. By fitting a NFW model to the reduced tangential shear profiles they concluded that the masses at $R_{500}^2$ are biased low by ∼7 − 10% with a scatter of 20%. The lensing analysis based on simulations in Becker & Kravtsov (2011) also concluded that, including shape noise only, these mass estimates are biased by ∼5 − 10% with a 20–30% scatter. In the same line, the weak lensing analysis in Euclid Collaboration (2024) derived an average 5% mass bias at R₂₀₀, with the bias for a given cluster differing by up to 30% depending on the orientation of the projection of the cluster observation.

Giocoli et al. (2012a) investigated also the dependence of the bias with redshift and found very little evolution. Quite the opposite, they showed that the radial range chosen to fit the density or mass profile model has an important impact on the bias. According to Euclid Collaboration (2024), the underestimation of cluster masses in Euclid-like weak lensing reconstructions is more important at higher redshift, modulated by the number density of background sources.

Some works compared also the uncertainties of masses reconstructed from total matter observables to those estimated from gas observables. They showed, on the one hand, the small impact of orientation in gas observables with simulated SZ and X-ray data (Buote & Humphrey 2012) and, on the other hand, the irreducible non-sphericity of haloes that affects the intrinsic scatter of the weak lensing masses, tracers of the total matter distribution (Becker & Kravtsov 2011). In the same line, Meneghetti et al. (2010b) suggested that lensing masses are three times more scattered (∼17 − 23% of scatter) than X-ray estimates. The error budgets given in Pratt et al. (2019) indicate that the assumption of spherical symmetry affects HSE masses at the level of a few percent, while lensing masses are affected about 10% due to mass modelling. According to Rasia et al. (2012), weak lensing biases are at least twice more scattered than X-ray outcomes.

In this work, we aim at quantifying and comparing the different contributions to the bias and its scatter. We perform the analysis in the framework of the NIKA2 Sunyaev-Zel’dovich Large Programme (LPSZ, Sect. 2.1) with simulated clusters from THE THREE HUNDRED project (Sect. 2.2.1, Cui et al. 2018). Given the importance of the sample selection, we use clusters that are twins to those from the LPSZ (Sect. 2.2.2), so that at least the mass and redshift distribution is representative of the real sample. We focus the work on total mass observables, but we also compare the results of gas observables of clusters. We compare all the masses at R₅₀₀ to be consistent with the LPSZ outcomes.

In a recent analysis in Euclid Collaboration (2024), the same THE THREE HUNDRED simulated data were used to study the bias of weak lensing mass reconstructions as part of the preparation work for the Euclid mission (Laureijs et al. 2011). Given the observational framework they wanted to mimic, in Euclid Collaboration (2024) the simulated data were converted into a Euclid-like weakly lensed galaxy distribution. With these galaxies, and accounting for the dispersion of the shapes of galaxies, they built the excess surface mass density profiles. The latter were fitted with the aforementioned tNFW model to reconstruct the three-dimensional mass density profiles of clusters. This paper does not aim to reproduce such observational effects and comparisons to Euclid Collaboration (2024) are valuable to elucidate the origin of some of the observed effects in that work. In addition, this paper differs from the work in Euclid Collaboration (2024) regarding the considered cluster sample and the overdensity of the reconstructed masses. While we restrict our analysis to a subsample of THE THREE HUNDRED clusters in the redshift range covered by the LPSZ programme and evaluate the masses at R₅₀₀, in Euclid Collaboration (2024) the authors studied the bias at R₂₀₀ for all THE THREE HUNDRED clusters at 9 different redshifts between z = 0.12 and z = 0.98. That being so, the complementarity to the analysis in Euclid Collaboration (2024) is clear.

This paper is organised as follows. In Sect. 2 we present the framework of this work: the NIKA2 LPSZ, THE THREE HUNDRED project and the clusters selection. The different data used are presented in Sect. 3. In Sect. 4 we describe the method followed to reconstruct the mass from convergence maps and the related bias. The different contributions to the scatter of the bias are discussed in Sect. 5. Finally, in Sect. 6 we compare the properties of mass biases obtained from total mass observables to gas observables used to reconstruct masses under hydrostatic equilibrium hypothesis. Conclusions are given in Sect. 7 together with the summary of our main results.

2. THE THREE HUNDRED-NIKA2 LPSZ samples

2.1. The NIKA2 Sunyaev-Zel’dovich Large Programme

The NIKA2 Sunyaev-Zel’dovich Large Programme or LPSZ (Perotto et al. 2022; Mayet et al. 2020) is a high angular resolution follow-up of ∼45 clusters of galaxies. These clusters were detected with Planck (Planck Collaboration XVIII 2015) or the Atacama Cosmology Telescope (ACT, Hasselfield et al. 2013) and were chosen to cover a wide range of masses (3 ≤ M₅₀₀/10M_⊙ ≤ 10) at intermediate to high redshift (0.5 ≤ z ≤ 0.9). The LPSZ benefits from high angular SZ observations with the NIKA2 camera (Adam et al. 2018; Bourrion et al. 2016; Calvo et al. 2016) and X-ray data obtained with the XMM-Newton and Chandra satellites.

By combining SZ and X-ray data, one of the objectives of the LPSZ is to re-estimate precisely the HSE M₅₀₀ mass of the galaxy clusters in the sample by reconstructing resolved hydrostatic mass profiles. This will allow us to build the relation between M₅₀₀ and the SZ signal for clusters at redshifts above 0.5 and to compare it to the relations measured at lower redshifts (Arnaud et al. 2010, hereafter A10).

However, the HSE hypothesis used to derive mass profiles is known to bias low the estimates (Pratt et al. 2019). This bias is currently being investigated with different approaches (Gianfagna et al. 2021; Salvati et al. 2019; de Andres et al. 2022), since increasing the value of the HSE mass bias is considered as a possible solution to release the tension between the cosmological results obtained from cluster counts and CMB measurements (Planck Collaboration XXI 2013; Planck Collaboration XXIV 2015; Salvati et al. 2018). To correct for such bias, mass measurements that do not rely on the hydrostatic equilibrium assumption are needed and lensing masses can be used for this purpose.

To tackle this issue for the NIKA2 LPSZ hydrostatic masses, we compared the HSE masses to lensing ones using the Cluster Lensing And Supernova survey with Hubble (CLASH, Postman et al. 2012) data for several clusters (Muñoz-Echeverría et al. 2022, 2023; Ferragamo et al. 2022). Our analyses were based on the publicly available convergence maps reconstructed from weak and strong lensing data (Zitrin et al. 2015). As described in Sect. 4.1, we fitted density models to those convergence maps, being able to get, from the best-fit density profiles, the lensing masses of clusters. In this work, we follow the same method on simulated clusters and we mainly focus our study on the biases that arise in the reconstruction of lensing masses from convergence maps.

2.2. Twin samples

2.2.1. THE THREE HUNDRED project

This work is based on THE THREE HUNDRED project⁴ galaxy cluster simulations. THE THREE HUNDRED project consists of hydrodynamical re-simulations of the 324 Lagrangian regions centred on the most massive galaxy clusters identified in the MultiDarkPlanck2 box of side length 1 h⁻¹ Gpc of the MultiDark⁵ dark-matter-only simulation. The regions, of radius 15 h⁻¹ Mpc and identified at z = 0, contain clusters with virial masses above 1.2 × 10¹⁵ M_⊙ and dark-matter particles of mass 1.5 × 10⁹ h⁻¹ M_⊙. The simulation assumes a cosmology based on the Planck Collaboration XIII (2016) results: h = 0.6777, n = 0.96, σ₈ = 0.8228, Ω_Λ = 0.692885, Ω_m = 0.30711, and Ω_b = 0.048206. This is also the cosmological model assumed in the rest of this paper.

The selected dark-matter-only volumes were mapped back to the initial conditions, dark matter particles were split into dark matter and gas particles and the regions were re-simulated with three hydrodynamical codes: GADGET-MUSIC (Sembolini et al. 2013), GIZMO-SIMBA (Davé et al. 2019; Cui et al. 2022), and GADGET-X (Rasia et al. 2015). We use the outputs coming from the latter, which contains very complete baryonic physics models (Table 2 in Cui et al. 2018) with dark matter, gas, stellar, and black hole particles. Previous works based on this code have shown the agreement between simulations and observations regarding: gas density and entropy profiles (Rasia et al. 2015), pressure profiles (Planelles et al. 2017), and gas density and temperature at around R₅₀₀ (Truong et al. 2018).

2.2.2. Twin samples definitions

For our analysis, we used the clusters from THE THREE HUNDRED project simulations selected to constitute the twin samples of the LPSZ. As described in Paliwal et al. (2022), the clusters were chosen to cover the same redshift range as the LPSZ, that is, 0.5 ≲ z ≲ 0.9. Amongst all the clusters in THE THREE HUNDRED satisfying this condition, we chose the snapshots 101 at z = 0.817, 104 at z = 0.700, 107 at z = 0.592, and 110 at z = 0.490. Three different samples were generated by matching properties of the clusters in the simulation to those known for the clusters in the LPSZ. The three twin samples are: (1) TS_{Tot − M}, the total mass twin sample, in which the clusters in the simulation are selected so that the total M₅₀₀ from the simulation matches with the M₅₀₀ mass of the LPSZ clusters according to Planck or ACT catalogues (Planck Collaboration XVIII 2015; Hilton et al. 2021), (2) TS_{HSE − M}, the hydrostatic mass twin sample, in which the hydrostatic mass of simulated clusters, $M_{500}^{\mathrm{HSE}}$, is matched with the Planck or ACT masses and (3) TS_Y, the twin sample based on the Y₅₀₀ parameter, the integrated SZ signal within R₅₀₀.

This selection gives three twin samples with 45 clusters (Fig. 2 in Paliwal et al. 2022). We present in Table 1 the range of masses covered by the clusters in each twin sample. Some of the clusters being repeated, we have altogether 122 different objects. In this work, we use all the clusters from the three twin samples simultaneously. The effects studied in the following do not vary from one twin sample to another and, on the contrary, accounting for the 122 clusters improves significantly our statistics.

3. The synthetic datasets

In this section we present the data generated from the simulations and exploited in our analyses.

3.1. Spherical thermodynamical profiles

For comparison sake, we make use of the spherical three-dimensional profiles of thermodynamical properties extracted from the simulations. In particular, we use the ICM pressure and total mass profiles. The profiles were computed accounting for particles in concentric shells and spheres (two-dimensional and three-dimensional profiles, respectively) centred at the maximum density peak, with radial bins starting at 10 kpc from the centre and increasing radii by 10% (as in Gianfagna et al. 2022).

The gas pressure profiles were extracted considering hot gas particles (i.e. only particles with temperature above 0.3 keV) and following Eq. (1) in Planelles et al. (2017), with the correction for SPH simulations proposed in Planelles et al. (2017) and Battaglia et al. (2015). The total mass profiles account for the mass of all the particles (dark matter, hot and cold gas, stars and black holes) and are used to define the true R₅₀₀ and M₅₀₀ of clusters. Profiles at r ≲ 0.025 × R₅₀₀ are not reliable since less than 100 particles are used to measure the thermodynamical quantities at these radial ranges.

In Fig. 1 we present the normalised pressure and mass density profiles for all the clusters in our sample. The pressure profiles have been normalised with respect to the P₅₀₀ obtained from Eq. (5) in A10. Following Fig. 3 in Gianfagna et al. (2021), we plot in Fig. 1 the “universal” pressure profile from A10 as well as the profile obtained from the combination of XMM-Newton and Planck data in Planck Collaboration Int. V (2013), hereafter P13. The density profiles were obtained by dividing the mass in each spherical shell by the volume of the shell and normalised with the critical density of the Universe at the corresponding redshift of each cluster². The panels in Fig. 1 show the self-similarity of the clusters regarding the gas pressure distribution, as well as the total matter density distribution. Clusters appear consistent with the self-similar hypothesis.

Fig. 1. Normalised gas pressure (left) and total matter mass density (right) spherical profiles for all the clusters in THE THREE HUNDRED-NIKA2 LPSZ twin samples. We also plot the “universal” pressure profile from the X-ray analysis in A10 and the X-ray and SZ analysis in P13. All profiles are drawn with respect to the radii normalised by the R500 corresponding to each cluster. The grey shaded areas indicate the regions below 0.025 × R500 that are not numerically robust.

3.2. Projected observable maps

3.2.1. Projected total mass: κ-maps

To study the projected total mass of galaxy clusters, we use the convergence maps or κ-maps produced as described in Herbonnet et al. (2022) and Euclid Collaboration (2024). These convergence maps were generated together with shear maps to mimic the lensing effect due to the simulated clusters on background sources. The convergence of a lens at a position θ is defined as

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })=\frac{\mathrm{\Sigma }(\mathit{\theta })}{{\mathrm{\Sigma }}_{\mathrm{crit}}},\end{array}$$(1)

where Σ(θ) is the projected mass density of the cluster at θ and Σ_crit the so-called critical surface density:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_S}{D_L{D}_{\mathit{LS}}}·\end{array}$$(2)

Here D_S, D_L, and D_LS are the angular diameter distances between the observer and the background source, the observer and the lens (the cluster), and the source and the lens, respectively.

To create the convergence maps, first the Σ-maps were obtained by projecting, within a volume of depth 10 Mpc, the masses of all the particles of the clusters along different axes. Then, these Σ-maps were divided by Σ_crit, assuming all background sources to be at z = 3. For each cluster, maps of 6  ×  6 Mpc² were produced, centred on the minimums of the gravitational potential wells. To avoid boundary errors, final maps are 5  ×  5 Mpc² with 2048  ×  2048 pixels, which gives different angular resolution pixels for different redshift clusters.

Six different maps are available for each cluster, corresponding to different projection axes. Three maps (hereafter 0, 1, and 2) are obtained from the projection along the x, y, and z orthogonal axes of the simulation. In principle, they correspond to random directions with respect to the morphology of the clusters. The other three (0_pr_axes, 1_pr_axes, and 2_pr_axes) are aligned with the principal axes of the clusters regarding their moments of inertia (Knebe et al. 2020).

3.2.2. Sunyaev-Zel’dovich effect: y-maps

The gas in the ICM of galaxy clusters is investigated starting from the y-maps. The thermal SZ effect due to galaxy clusters is characterised by the Compton parameter, proportional to the integral along the line of sight of the thermal pressure P_e in the ICM,

$$\begin{array}{c}\hfill y=\frac{{\sigma }_{\mathrm{T}}}{m_e{c}^2}{\displaystyle {\int }_0^{+\infty }}{P}_e\mathrm{d}l=\frac{{\sigma }_{\mathrm{T}}{k}_{\mathrm{B}}}{m_e{c}^2}{\displaystyle {\int }_0^{+\infty }}{n}_e{T}_e\mathrm{d}l,\end{array}$$(3)

with σ_T, k_B, m_e, and c the Thomson cross section, the Boltzmann constant, the electron rest mass, and the speed of light, respectively. Assuming the ICM behaves as an ideal gas, the Compton parameter is also directly proportional to the integral of the electron number density n_e and the electron temperature T_e.

As described in Cui et al. (2018) and previously done in Sembolini et al. (2013) and Le Brun et al. (2015), in THE THREE HUNDRED project this integration is converted into a summation⁶ over the gas particles in the line of sight and within a 2 × R₂₀₀ depth volume,

$$\begin{array}{c}\hfill y=\frac{{\sigma }_{\mathrm{T}}{k}_{\mathrm{B}}}{m_e{c}^2\mathrm{d}A}{\displaystyle \sum _{i=0}^{N_{\mathrm{gas}\mathrm{part}.}}}{T}_{e,i}{N}_{e,i}W(r,h_i).\end{array}$$(4)

Here the electron number density n_e is represented as the number of electrons in a given gas particle, N_e, divided by its spatial volume dV (dV = dldA): n_e = N_e/dV = N_e/dA/dl. In Eq. (4) W(r, h_i) corresponds to the SPH smoothing kernel used to smear the signal of each particle to the projected pixels, with smoothing length h_i (Cui et al. 2018).

For each cluster, there are 29 y-maps available, projected along the x, y, and z main axes of the simulation, plus along other 26 random axes. We analysed only the projections in common with the κ-maps, hence the three main axes projections. Centred on the projected maximum density peak of the clusters, maps have 1920 × 1920 pixels of 5′′ angular resolution, going up to 14 − 30 × R₂₀₀.

3.2.3. Projected gas mass: Gas mass maps

A different way to characterise the distribution of the gas in the ICM is to use projected gas mass maps. They were generated as the Σ-maps (Sect. 3.2.1), but accounting only for the gas particles along the line of sight. Final maps are also 5 × 5 Mpc² with 2048 × 2048 pixels and projected along the three main axes of the simulation.

4. Total mass reconstruction from κ-maps

Just as in Meneghetti et al. (2014), in this work we do not consider any systematic effect, such as the accuracy of galaxy shape measurements or photometric redshifts, that may affect lensing analyses. Those effects can introduce additional uncertainties that require separate investigations (Bahé et al. 2012; Euclid Collaboration 2024). We restrict our work to the effects that arise from the reconstruction of three-dimensional quantities from projected data.

4.1. Mass reconstruction procedure

As aforementioned, to estimate the total mass of the clusters from projected maps we followed the procedure used in Ferragamo et al. (2022) and Muñoz-Echeverría et al. (2022, 2023). In these works convergence maps from the CLASH dataset are used to reconstruct the lensing mass of galaxy clusters: κ-maps are converted into projected mass density maps, Σ, and mass density models are fitted to the radial profile of Σ, which allows one to reconstruct the lensing mass.

From the convergence maps described in Sect. 3.2.1, we obtain projected Σ-profiles by radially averaging the values in the map, starting from their central pixel. Uncertainties are computed from the dispersion in each radial bin to account for non-circular features in the map. In a previous work Giocoli et al. (2012a) showed the important impact of the chosen radial range to perform the fit. To build the excess surface mass density profiles in Euclid Collaboration (2024) the authors used 22 logarithmically spaced bins, spanning from 0.02 to 1.7 h⁻¹ Mpc from the cluster centre.

We decided to consider instead a radial range in which the fractional errors of the Σ-profiles are comparable to the uncertainties obtained from data-based convergence maps. In the left panel in Fig. 2, we show the variance of binned profiles obtained from convergence maps. In grey we present the profiles obtained from the simulated convergence maps. We aim to assess the lensing mass reconstruction of the NIKA2 LPSZ clusters (Ferragamo et al. 2022; Muñoz-Echeverría et al. 2022, 2023) based on the CLASH convergence maps (Zitrin et al. 2015) obtained with the highest quality lensing data. Thus, we show in red, green, blue, and yellow in Fig. 2 the variance of the convergence profiles obtained from CLASH convergence maps for the four clusters (with two convergence map models for three of them) analysed in Muñoz-Echeverría et al. (2022, 2023) and Ferragamo et al. (2022; for more details on the uncertainties of the profiles see Muñoz-Echeverría et al. 2023). We observe an overall agreement, with possible differences in the variance trends in regions < 0.02 arcmin, where all simulation profiles tend to rapidly decrease.

Fig. 2. Variance of the convergence profiles (left) and fractional error of the Σ-profiles (right) calculated from convergence maps. In grey we show the profiles for the simulated clusters in THE THREE HUNDRED-NIKA2 LPSZ samples. Red, green, blue, and yellow profiles correspond to CL J1226.9+3332, PSZ2 G144.83+25.11, PSZ2 G228.16+75.20, and MACS J1423.8+2404 clusters, respectively, with data from the CLASH convergence maps (Zitrin et al. 2015) and uncertainties computed following Muñoz-Echeverría et al. (2023; sample defined in Muñoz-Echeverría et al. 2022). Solid and dashed lines differentiate the profiles for the CLASH convergence maps reconstructed assuming the Light-Traces-Mass (LTM) and the Pseudo-Isothermal Elliptical Mass Distribution with an elliptical NFW (PIEMD+eNFW) models, respectively. Vertical dotted lines in the right panel indicate the minimum and maximum of the used radial range.

In the right panel in Fig. 2 we compare the fractional errors of the Σ-profiles from the simulation and from CLASH maps. The agreement is again very good, but in the region between ∼0.02 and ∼0.05 arcmin the observed errors tend to be larger than the average simulated errors. For this reason, we restricted the fit of the Σ-profiles to the radii above 0.1′, removing the radial ranges where the fractional uncertainties of simulated Σ-profiles are systematically smaller than data-based error bars and, therefore, unrealistic. Regarding the upper limit, we restricted the fit of the profiles to 2′, which is the maximum radial extent of convergence maps from Zitrin et al. (2015).

To model the Σ-profile we used the two-parameter NFW mass density model:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{NFW}}(r)=\frac{{\rho }_s}{(r/r_s){(1+r/r_s)}^2}=\frac{c_{\mathrm{\Delta }}^3\mathrm{\Delta }{\rho }_{\mathrm{crit}}/3I(c_{\mathrm{\Delta }})}{(r{c}_{\mathrm{\Delta }}/R_{\mathrm{\Delta }}){(1+r{c}_{\mathrm{\Delta }}/R_{\mathrm{\Delta }})}^2},\end{array}$$(5)

and the three-parameter Einasto model,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{Ein}}(r)={\rho }_s\exp \{-\frac{2}{\alpha }[{\left(\frac{r}{r_s}\right)}^{\alpha }-1]\},\end{array}$$(6)

and fitted their projected mass density profiles to the Σ-profiles from the convergence maps:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{model}}(R)=2{\displaystyle {\int }_0^{+\infty }}{\rho }_{\mathrm{model}}\left(\sqrt{r^2+R^2}\right)\mathrm{d}r.\end{array}$$(7)

Here r is the radius of three-dimensional spherical profiles while R is used for the projected profiles. In Eq. (5), c_Δ is the concentration parameter that corresponds to R_Δ, c_Δ = R_Δ/r_s, and I(c_Δ) is a function that depends only on c_Δ and the shape of the density profile (Eq. (15) in Aguena et al. 2021), which for the NFW model gives

$$\begin{array}{c}\hfill I(c_{\mathrm{\Delta }})=ln(1+c_{\mathrm{\Delta }})-\frac{c_{\mathrm{\Delta }}}{1+c_{\mathrm{\Delta }}}·\end{array}$$(8)

The parameters for the Einasto density model (Eq. (6)) are r_s, ρ_s, and α: the scale radius, the characteristic density, and the shape parameter, respectively.

The fits were performed via a Markov chain Monte Carlo (MCMC) analysis using the emcee Python package (Foreman-Mackey et al. 2019; Goodman & Weare 2010). In addition, we made use of the profiley Python package⁷, which contains already tested (Madhavacheril et al. 2020) functions that describe NFW and Einasto spherical density profiles, as well as the corresponding line of sight projected profiles. Regarding the prior ranges of the model parameters, we considered uniform priors (as in Ettori et al. 2019) wide enough to cover all possible results found in the literature. For the NFW model we took c₂₀₀ = 𝒰(0.1, 15) and r_s = 𝒰(0.1, 6) Mpc, which are safe ranges in view of the concentration and mass values obtained for observational results in Ettori et al. (2019) and Umetsu et al. (2020). For Einasto we considered r_s = 𝒰(0.05, 3) Mpc, ${\rho }_s=\mathcal{U}({10}^{10},{10}^{18})M_{\odot }{\mathrm{Mpc}}^{-3}$, and α = 𝒰(0, 5) to make sure to cover broadly the variety of possible results (Ettori et al. 2019; Eckert et al. 2022b). Once the fit was performed, we verified that the modelled Σ-profile is a good representation of the data following a goodness of fit criterion defined as

$$\begin{array}{c}\hfill \mathrm{Fit}\mathrm{quality}=\mathrm{median}(|{\mathrm{\Sigma }}_{\mathrm{data}}-{\mathrm{\Sigma }}_{\mathrm{model}}|/{\sigma }_{{\mathrm{\Sigma }}_{\mathrm{data}}}),\end{array}$$(9)

where the median is estimated over all the radial bins and posterior distributions of the model. Bad fits were rejected if Fit quality >  1. This means that if the median of the absolute difference between the Σ-profile data bins and the model with respect to the uncertainties is larger than unity, the fit is not satisfactory.

After the fit, mass profiles were obtained by computing for each sample of parameters the corresponding profile. For NFW model, the spherical mass profile is described by,

$$\begin{array}{cc}\hfill M_{\mathrm{NFW}}(<r)& =4\pi {\rho }_s{r}_s^3[\frac{1}{1+r/r_s}+ln(1+r/r_s)-1]\hfill \\ \hfill & =4\pi c_{200}^{3}200{\rho }_{\mathrm{crit}}/3I(c_{200}){(R_{200}/c_{200})}^3\hfill \\ \hfill & [\frac{1}{1+r{c}_{200}/R_{200}}+ln(1+r{c}_{200}/R_{200})-1].\hfill \end{array}$$(10)

And for Einasto,

$$\begin{array}{c}\hfill M_{\mathrm{Ein}}(<r)=\frac{4\pi {\rho }_s{r}_s^3}{\alpha }{e}^{2/\alpha }{\left(\frac{\alpha }{2}\right)}^{3/\alpha }\gamma [\frac{3}{\alpha },\frac{2}{\alpha }{\left(\frac{r}{r_s}\right)}^{\alpha }],\end{array}$$(11)

where γ is the incomplete lower gamma function:

$$\begin{array}{c}\hfill \gamma (a,x)={\displaystyle {\int }_0^x}{t}^{a-1}{e}^{-t}\mathrm{d}t.\end{array}$$(12)

Finally, from the reconstructed mass profiles we computed the $M_{500}^{\kappa }$ mass of the cluster at its corresponding $R_{500}^{\kappa }$ (Fig. A.2). This procedure was repeated for each κ-map, therefore, several times for each cluster, which gave different $M_{500}^{\kappa }$ estimates affected by the projection effect and the choice of the density model. As an example, in Fig. 3 we present the Σ-profiles for the six projections of one cluster (number 0306) at z = 0.817. The orange and purple lines show the mean profiles (solid lines) and 1σ error bars (shaded areas) obtained from the posterior distributions of NFW and Einasto fit parameters. In the 2_pr_axes projection the cluster is observed along its most elongated axis. For this reason, the centre appears very massive and the density in the outskirts drops fast. We verified that the sizes of the uncertainties for these Σ-profiles are compatible with the ones obtained from the convergence maps built by combining weak and strong lensing data. Nevertheless, they do not necessarily correspond to what we would expect from other observations. For instance, in weak lensing mass density reconstructions error bars tend to be larger in the centre and smaller in the outskirts of clusters.

Fig. 3. Projected mass density profiles (empty markers) for the six convergence maps for the 0306 cluster in the snapshot 101 (z = 0.817). The orange and purple profiles show the fit with NFW and Einasto models, respectively. We give the mean profiles with 1σ contours. From left to right and top to bottom we present the 0, 1, and 2 random projections and 0_pr_axes, 1_pr_axes, and 2_pr_axes principal moment of intertia projections. The vertical cyan line shows the true R500 of the cluster.

We repeated the analysis for all the clusters in our sample and for the following we take, amongst all the Σ-profile fits, the results that pass the quality criteria in Eq. (9) for both the NFW and Einasto fits. As explained in Sect. 5.2, in the rest of the study we only consider the three projections along the x, y, and z orthogonal axes of the simulation.

We study the correlations between the best-fit values of the posterior distributions of model parameters for all the clusters (shown in Fig. A.1). An anticorrelation between Einasto’s ρ_s and r_s parameters is clearly noticeable, with a Pearson correlation coefficient of −0.56. Contrary to some studies (López-Cano et al. 2022, and references therein) showing a decreasing tendency of the concentration, c₂₀₀, with the mass, M₂₀₀, we do not find any significant correlation. In Cui et al. (2018) the concentration-mass relation of z = 0 clusters from THE THREE HUNDRED simulation is presented, showing that the relation is flatter for GADGET-X clusters than for GADGET-MUSIC. Some works (e.g. Oguri et al. 2009) also show a dependence of the concentration parameter with redshift, but our results do not present any clear evolution along cosmic time (the Pearson correlation coefficient between the redshift and the concentration is of r = −0.09). According to Meneghetti et al. (2014), the redshift evolution of the concentration-mass relation is weak for the MUSIC-2 simulation.

4.2. Mass bias from convergence maps

We define the mass bias as the relative distance of each reconstructed mass, $M_{500}^{\kappa }$ (as above-mentioned, evaluated at the $R_{500}^{\kappa }$ obtained from the reconstructed mass profile), to the true mass of the cluster, M₅₀₀:

$$\begin{array}{c}\hfill b^{\kappa }=(M_{500}-M_{500}^{\kappa })/M_{500}.\end{array}$$(13)

We present in Fig. 4 the bias estimated from all the analysed maps with respect to the true M₅₀₀ for each cluster. Orange and purple markers show the bias and uncertainties measured for the NFW and Einasto mass density models, respectively, with the uncertainties from the fits. The horizontal orange line shows the mean bias for NFW, while the purple corresponds to Einasto. Einasto gives, on average, a less biased estimate of the mass of the clusters. The median biases for NFW and Einasto reconstructions are med(b^κ)=0.097 and 0.086, respectively. We identify the origin of the bias from a thorough examination of the Σ-profiles presented in Sect. 4.1.

Fig. 4. Mass bias for masses inferred from fitted κ-maps with respect to M500 per cluster. Each data point corresponds to the mean bias of a cluster for a given random projection (0, 1 or 2) with error bars showing the 1σ uncertainty. Results for the Einasto model have been artificially shifted in mass for visualisation purposes. The horizontal orange and purple dashed lines show the mean for all the results for NFW and Einasto, respectively. The black line indicates the zero.

We present in Fig. 5 the relative difference between the fit of the Σ-profile with NFW and Einasto models and the profiles measured from the convergence maps. We observe that for both models the mass density profiles are systematically underestimated for the radial ranges close to R₅₀₀. This is then translated into biased mass profiles and, consequently, biased M₅₀₀ estimates.

Fig. 5. Relative differences of Σ-profiles between the fits with NFW (top) and Einasto (bottom) models and the profiles measured from the convergence maps. Solid lines show the median.

One could think of improving the fits either by adapting the radial range ad hoc to get a less biased result or by choosing a more flexible density model. As aforementioned, this type of studies have already been done in the literature (Meneghetti et al. 2014; Giocoli et al. 2012a) and we do not go through those refinements in this work. In Cui et al. (2018), the authors concluded that the NFW model is a good representation of the spherical mass density profile of the clusters in THE THREE HUNDRED simulation, even considering different inner radii. But in Euclid Collaboration (2024), the authors showed that THE THREE HUNDRED clusters do not follow perfect NFW density profiles.

In the left panel in Fig. 6 we present the mean bias per cluster obtained from the NFW density fits with respect to the mean bias from the Einasto fits. There is a strong correlation between the NFW and Einasto biases, with a Pearson correlation coefficient of 0.86. This confirms that the impact of the model choice is not enough to blur the information in the convergence maps.

Fig. 6. Mean bias of the three projections (left) and ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ (right) per cluster for Einasto against NFW reconstructions. Each dot corresponds to one cluster. We give the Pearson correlation coefficient for the values given in each panel.

We present in Fig. 7 the bias measured from the maps projected along the three axes of inertia and verify that the orientation of the cluster affects directly the bias: when the cluster is elongated along the line of sight (2_pr_axes) the reconstructed mass is on average overestimated and when the major axes are on the plane of the sky (0_pr_axes), underestimated. In Herbonnet et al. (2022), the authors investigated the relation between the orientation of haloes and their brightest cluster galaxy (BCG). They found that the major axes of haloes and BCGs are aligned with an average separation of ∼20°. For this reason, having access to the orientation of the BCG would provide a way to improve the knowledge of the orientation of the cluster, and therefore, of its mass reconstruction.

Fig. 7. Mass bias for NFW (orange) and Einasto (purple) fits of the Σ-profiles when aligned with the principal axes of the clusters regarding their moments of inertia. Solid lines correspond to 0_pr_axes, dashed lines to 1_pr_axes, and dotted lines to 2_pr_axes. We give the value of the median and the 16th and 84th percentiles. Shaded areas represent the 1σ regions of the bias obtained for all the clusters with the three random projections (Fig. 4).

Given all the mentioned differences between our work and the analysis in Euclid Collaboration (2024), comparisons to the latter have to be done with extraordinary care. By fixing the truncation parameter of the tNFW model to t = 3, they obtained that the weak lensing masses reconstructed for the clusters at redshift z = 0.22 are on average biased low by ∼7%, with a standard deviation of 0.24 over the full sample (using the full THE THREE HUNDRED sample of clusters). Assuming also a fixed concentration parameter c₂₀₀ = 3 and considering the clusters at redshifts below z = 0.7, according to Euclid Collaboration (2024), weak lensing M₂₀₀ masses are on average underestimated by less than 5%. We observe in Fig. 14 in Euclid Collaboration (2024) that, for the redshift range considered in our work (0.5 <  z <  0.9), the mean biases on M₂₀₀ vary from 0% to 20% depending on the chosen value for the c₂₀₀ parameter. These results match our findings at R₅₀₀. Regarding the impact of the orientation of clusters on the mass bias, in Euclid Collaboration (2024) the authors concluded that when clusters are oriented along the major and minor axes, the mass biases are ∼25% larger and smaller, respectively, with respect to the biases for random orientations. Therefore, their results are completely in line with Fig. 7.

5. Scatter of the mass bias

In the previous section we have been looking at the mean biases estimated with each model for our full sample. The dispersion from cluster to cluster and for the different projections of the same cluster appear to be very important (Fig. 4). In this section we try to identify and quantify the different effects that contribute to the scatter of the biases in the reconstruction. For this, we use an statistical approach and we consider the full sample.

We first define the total variance per model from the scatter of the bias values obtained with the same model for all the clusters and all their projections (see Fig. 4):

$$\begin{array}{c}\hfill {\sigma }_{b,\mathrm{tot}}^2=\frac{1}{N_{\mathrm{chains}}{N}_{\mathrm{projections}}{N}_{\mathrm{clusters}}}{\displaystyle \sum _{i,j,k=1,1,1}^{N_{\mathrm{chains}},N_{\mathrm{projections}},N_{\mathrm{clusters}}}}{(b_{i,j,k}^{\kappa }-⟨b^{\kappa }⟩)}^2,\end{array}$$(14)

where N_chains is the amount of chains kept after each fit, N_projections is the number of different projections per cluster and N_clusters the amount of clusters considered. Here ⟨b^κ⟩ is the mean bias for all the clusters, accounting for all the projections and chains.

We separate this total variance in two main contributions: (1) the variance from cluster to cluster, ${\sigma }_{b,\text{cluster-to-cluster}}^2$, showing how different are the bias values depending on the cluster, and (2) the variance of the results for each cluster, due to the different values depending on the projection, ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ . So we write (as in Bartalucci et al. 2023):

$$\begin{array}{c}\hfill {\sigma }_{b,\mathrm{tot}}^2\sim {\sigma }_{b,\mathrm{cluster}-\mathrm{to}-\mathrm{cluster}}^2+{\sigma }_{b,3\mathrm{projections}}^2.\end{array}$$(15)

Here ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ is the dispersion of the results of three projections, that for each cluster c is obtained as,

$$\begin{array}{c}\hfill {({\sigma }_{b,3\mathrm{projections}}^2)}_c=\frac{1}{N_{\mathrm{chains}}{N}_{\mathrm{projections}}}{\displaystyle \sum _{i,j=1,1}^{N_{\mathrm{chains}},N_{\mathrm{projections}}}}{(b_{i,j,c}^{\kappa }-⟨b^{\kappa }{⟩}_c)}^2,\end{array}$$(16)

where ⟨b^κ⟩_c is the mean bias per cluster c, considering the different projections and all the chains.

At the same time, in each ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ per cluster there are two sources of dispersion that we name: (1) ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$, which accounts for the fact that the results obtained from three projections differ, and, (2) ${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2$ that quantifies the uncertainties of the fits. Therefore, we can write

$$\begin{array}{c}\hfill {\sigma }_{b,3\mathrm{projections}}^2\sim {\sigma }_{b,\mathrm{intrinsic}\mathrm{proj}}^2+{\sigma }_{b,\mathrm{fit}}^2.\end{array}$$(17)

In the following subsections we study each contribution to the total scatter.

5.1. Fit uncertainties: ${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2$

The most evident source of scatter is the one related to the model fitting uncertainties, that is, the scatter propagated from the posterior distributions of the parameters in the fit of the convergence map at one specific projection of one cluster. We define (${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2$)_proj, c as the variance of the bias values of each map fit (“proj”) for each cluster (c):

$$\begin{array}{c}\hfill {({\sigma }_{b,\mathrm{fit}}^2)}_{\mathrm{proj},c}=\frac{1}{N_{\mathrm{chains}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{chains}}}}{(b_{i,\mathrm{proj},c}^{\kappa }-⟨b^{\kappa }{⟩}_{\mathrm{proj},c})}^2.\end{array}$$(18)

We give the mean and median values in Table 2 (see Fig. B.1 for their distribution). On average 8% (${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2\sim 0.007$) and 9% (${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2\sim 0.008$) error in the mass reconstruction comes only from the uncertainty in the fit for the NFW and Einasto models, respectively.

5.2. Scatter from three projections: ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$

For each cluster three different projections are available and the variance of the results from the three is the ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ defined in Eq. (16). As for the fit uncertainties, we give the mean and median values in Table 2 (the histograms of the variance for each cluster, accounting for the posterior distributions of the three projections together are shown in Fig. B.1). In the right plot of Fig. 6 we present the ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ obtained with NFW and Einasto for each cluster in our sample, showing the correlation between the variances measured with the different models.

We study this scatter to quantify the uncertainty that should be added to observational mass reconstructions for which clusters are only observed in one projection. To truly estimate the projection effect, the dispersion of the masses obtained from the infinite projections of each cluster would be needed. In Appendix C we instead compared the scatter obtained from 100 projections with the scatter obtained from the 3 simulation axes and from the 3 inertial axes and found that using the projections along x, y, and z (i.e. 0, 1, and 2) returns a more similar scatter level to the 100 projections. The projections along the inertia axes instead exhibit a significant larger scatter. For this reason, in the scatter calculations we did not make use of the maps projected along the inertia axes.

5.3. Intrinsic projection effect: ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$

For both NFW and Einasto models, there is a contribution to ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ that is not explained by the uncertainties on the fits. We consider that this is due to the real “intrinsic projection effect” presented in Eq. (17). From the difference between the mean (and also between the median) values of ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ and ${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2$, the intrinsic projection effect is of the order of ∼0.020 − 0.025. This corresponds to about 14 − 16% error on the mass bias and then on M₅₀₀. We stress that this is true for both NFW and Einasto, which is comforting, since we are looking for an intrinsic effect that should depend only on the clusters themselves.

Model independent intrinsic projection effect. We are also interested in quantifying the contribution of ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$ without being affected by modelling effects. For this purpose, we take a different approach and infer masses directly from the projected convergence maps. By integrating the Σ-maps up to the true θ₅₀₀, we compute directly the total mass of clusters from cylindrical integration, $M_{500}^{\kappa ,\mathrm{cyl}}$, where tan(θ₅₀₀)=R₅₀₀/𝒟_A, with R₅₀₀ the true radius of the cluster and 𝒟_A the angular diameter distance at the cluster redshift. For each cluster we have three different $M_{500}^{\kappa ,\mathrm{cyl}}$, one per projected map. The grey histogram in Fig. 8 shows the variance of the three mass biases per cluster c for all the clusters in our sample obtained from

$$\begin{array}{c}\hfill {({\sigma }_{b,\mathrm{intrinsic}\mathrm{proj}}^{\kappa -\mathrm{map}})}_c^2=\frac{1}{N_{\mathrm{projections}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{projections}}}}{(b_{j,c}^{\kappa -\mathrm{map}}-⟨b^{\kappa -\mathrm{map}}{⟩}_c)}^2.\end{array}$$(19)

Fig. 8. Variance of the mass bias estimated from maps. The grey histogram is obtained from the integration of convergence maps, i.e. $M_{500}^{\kappa ,\mathrm{cyl}}$. The magenta histogram corresponds to the variance for gas masses ($M_{500}^{\mathrm{Gas},\mathrm{cyl}}$) estimated from gas mass maps. The empty green histogram shows the dispersion of the $M_{500}^{\mathrm{cyl}}$ after conversion from the $Y_{500}^{\mathrm{cyl}}$ measured from y-maps, using the scaling relation in Fig. E.1.

On average, the projection effect when integrating the total mass at a fixed aperture is of the order of ${({\sigma }_{b,\mathrm{intrinsic}\mathrm{proj}}^{\kappa -\mathrm{map}})}^2\sim 0.01$, therefore, an error of 10% on the mass. The projection effect computed in this way is close, but not enough, to fully explain the ∼15% uncertainties obtained from the ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ − ${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2$ difference.

An additional term could explain the extra dispersion: while the masses used to compute ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ are estimated at different radii (each $M_{500}^{\kappa }$ measured at the corresponding $R_{500}^{\kappa }$), here all the $M_{500}^{\kappa ,\mathrm{cyl}}$ for ${({\sigma }_{b,\mathrm{intrinsic}\mathrm{proj}}^{\kappa -\mathrm{map}})}^2$ are measured at fixed θ₅₀₀ per cluster. Therefore ${({\sigma }_{b,\mathrm{intrinsic}\mathrm{proj}}^{\kappa -\mathrm{map}})}^2$ needs a corrective term to be fully comparable to ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$. To quantify this effect, we produced mock mass density profiles of galaxy clusters with a large variety of shapes and masses. We compared the scatter when estimating the M₅₀₀ spherically and up to the R₅₀₀ corresponding to each profile and the scatter when measuring the mass cylindrically integrated at a fixed aperture. All the details are given in Appendix D. In conclusion, this additional variance can be estimated of the order of ∼0.001 − 0.003, so ∼3 − 5%, depending on the shape, redshift, and mass of the density profile of the cluster. We note that additional differences could appear when integrating spherically against cylindrically, related, for example, to the presence of structures.

Therefore, both from a model dependent and a model independent approach we obtain compatible values for ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$. We conclude that when deriving cluster masses from observables tracing the total mass (for instance, the lensing effect), the intrinsic effect of the projection creates an error in the M₅₀₀ estimates between 10% and 16%. These values could vary with the depth of the considered κ-maps.

5.4. Total scatter: ${\sigma }_{b,\mathrm{t}\mathrm{o}\mathrm{t}}^2$

Finally, we quantify the total variance of the bias along the sample and the scatter from cluster to cluster. We compute, following the Eq. (14), the total variance of the bias in our sample by accounting for all the clusters, with all the projections and their fit uncertainties. Following Eq. (15), we obtain the excess with respect to the projection and fit effects. All the values are given in Table 2.

The total dispersion appears to be larger than the average projection effect. This could be due to some clusters with a very important projection effect (right panel in Fig. 6) that might be affecting significantly ${\sigma }_{b,\mathrm{t}\mathrm{o}\mathrm{t}}^2$. But it could also originate from the scatter of the mean bias values per cluster. We compute the mean bias per cluster c accounting for the masses reconstructed from all the chains of the three convergence maps ($⟨b^{\kappa }{⟩}_c=<(M_{500}-M_{500}^{\kappa })/M_{500}>$) and quantify how different are the mean biases for all the clusters in the sample. We define the variance of the mean biases, ${\sigma }_{b,\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}^2$, as:

$$\begin{array}{c}\hfill {\sigma }_{b,\mathrm{mean}}^2=\frac{1}{N_{\mathrm{clusters}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{clusters}}}}{(⟨b^{\kappa }{⟩}_k-⟨⟨b^{\kappa }{⟩}_c⟩)}^2,\end{array}$$(20)

where ⟨⟨b^κ⟩_c⟩ is the mean over the mean biases of all the clusters. Results are shown in Table 2. The missing dispersion in ${\sigma }_{b,\mathrm{t}\mathrm{o}\mathrm{t}}^2$ seems to be explained by the different mean bias values of different clusters, ${\sigma }_{b,\text{cluster-to-cluster}}^2$ ∼ ${\sigma }_{b,\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}^2$. We notice that it is a model dependent effect, NFW results being slightly more scattered with respect to ⟨⟨b^κ⟩_c⟩. This difference could mean that NFW reconstructions fail more often in the description of the mass of clusters and, compared to Einasto, NFW fit uncertainties are not large enough to account for it.

5.5. Intrinsic scatter: ${\sigma }_{b,\text{cluster-to-cluster}}^2$

The dispersion of the bias from cluster to cluster can be caused by multiple effects. We investigate here the impact of data limitations and of the mass, redshift, and dynamical state of the clusters.

Data limitations. We only consider 3 projections per cluster and the mean bias obtained from these 3 may not be a good representation of the true mean bias. We test the impact of this effect using four clusters for which we have 100 projections (see Appendix C). For each cluster we calculate the mean bias accounting for the 100 projections. Values are shown as vertical dashed lines in Fig. 9. At the same time, we take randomly three projections out of the 100 and compute the mean for those trios. The histograms in Fig. 9 show the results for 100 realisations. They illustrate that the mean bias ⟨b^κ⟩_c estimated from 3 random projections can in some cases be significantly different from the true one. Although histograms are centred in the mean bias of 100 projections, their variance is of the order of ∼0.005–0.010, which also contribute to ${\sigma }_{b,\text{cluster-to-cluster}}^2$ or ${\sigma }_{b,\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}^2$.

Fig. 9. Mean bias per cluster. The vertical lines show the mean bias when using 100 projections to reconstruct the mass. The histograms are obtained from the mean bias computed by taking randomly 3 projections out of the 100 available 100 times.

Fig. 10. Bias for clusters in different redshift (left) and mass (right) bins. We give the mean value and the σb, mean scatter for the mean bias of clusters in each bin. Results for NFW and Einasto models are shown in orange and purple, respectively, with a line as a guide to the eye.

Impact of mass and redshift. There could be also a correlation between intrinsic properties of clusters (mass, redshift, etc.) and the measured biases, meaning different responses to the model fitting depending on cluster properties. To investigate these effects, we separate our sample in subsamples. We choose a binning close to the one used for the NIKA2 LPSZ (Sect. 2.1): four bins in redshift (one per snapshot) and five bins in mass (in the NIKA2 LPSZ only two bins in redshift are considered, but with the objective of observing any evolution, we choose to distinguish four). The five bins in mass are: M₅₀₀/10¹⁴ M_⊙ < 4, 4 ≤ M₅₀₀/10¹⁴ M_⊙ < 6, 6 ≤ M₅₀₀/10¹⁴ M_⊙ < 8, 8 ≤ M₅₀₀/10¹⁴ M_⊙ < 10, and 10 ≤ M₅₀₀/10¹⁴ M_⊙.

We show in Fig. 10 the evolution of the mean bias with redshift and mass for NFW and Einasto mass reconstructions. A slight evolution of the bias with redshift is observed for NFW results, but not for Einasto. We note, in any case, that our results are compatible with no evolution. In Euclid Collaboration (2024), the authors observe an evolution of the bias of M₂₀₀ masses with redshift and claim that such evolution is due to the low weak lensing signal-to-noise at high redshift. In the right panel in Fig. 10, we observe no evolution of the bias with the true M₅₀₀ either.

Impact of the dynamical state. The dynamical state of clusters can also be an indicator of the bias of masses reconstructed from projected maps and their dispersion. To classify the clusters in our sample according to their dynamical state, we use the f_s and Δ_r indicators at R₂₀₀ (and R₅₀₀) that correspond respectively to the ratio between the sub-halo masses and the cluster mass within R₂₀₀ (R₅₀₀) aperture and the offset between the centre of mass of the cluster and the maximum density peak position, normalised by R₂₀₀ (R₅₀₀). Following the definitions in De Luca et al. (2021), clusters are relaxed (disturbed) if f_s <  0.1 and Δ_r <  0.1 (f_s >  0.1 and Δ_r >  0.1). Otherwise, we consider them in an intermediate dynamical state.

In the top panels in Fig. 11 we present the distribution of the mean bias per cluster independently for relaxed (green), intermediate (blue), and disturbed (red) clusters, classified with the indicators at R₂₀₀. The variance of these mean biases from cluster to cluster is larger for the disturbed sample (${\sigma }_{b,\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}^2$ = 0.0176 and 0.0070 for NFW and Einasto, respectively) than for the relaxed one (${\sigma }_{b,\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}^2$ = 0.0048 for NFW and 0.0055 for Einasto), but, on average, relaxed clusters are not less biased than disturbed ones.

Fig. 11. Bias properties depending on the dynamical state of clusters within R200. Top: distributions of the mean bias per cluster for relaxed (green), intermediate (blue), and disturbed (red) clusters. We show the bias from NFW (left) and Einasto (right) reconstructions. Bottom: distributions of the variance of the bias per cluster accounting for three projections. Clusters are classified according to their dynamical state and results for NFW (left) and Einasto (right) are shown. The vertical lines show the mean values of the distributions.

In the bottom panels in Fig. 11 we show the variance of the bias per cluster as computed from the three projections. Although we would expect to detect the departure from sphericity in the measurement of ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$, from these results we only observe such effect for the Einasto model. The three subsamples give scatter values that are compatible with those of the full sample.

The classification of clusters changes if we use the indicators at R₅₀₀. In Fig. 12 we summarise the mean and median values obtained from the distributions in Fig. 11 for dynamical states of clusters defined within R₅₀₀ and R₂₀₀. From this figure, the bias of masses reconstructed with our method seems vaguely related to the relaxation status of the galaxy clusters.

Fig. 12. Summarised bias properties depending on the dynamical state of clusters within R200 (empty) and R500 (filled). Left: mean (circle) and median (star) of the mean biases of each type of cluster. Right: mean and median of the ${\sigma }_{b,3\mathrm{projections}}^2$ for all the clusters of each type. We show in green, blue, and red the results for relaxed, intermediate, and disturbed clusters, respectively.

In Euclid Collaboration (2024), the authors used a stronger constraint to classify clusters as relaxed: they considered relaxed systems if f_s <  0.1 and Δ_r <  0.05 (x_off in Euclid Collaboration 2024) or if f_s <  0.1, Δ_r <  0.04, and 0.85 <  η <  1.15 (with η the virial ratio, as in Cui et al. 2018), all within R₂₀₀, and unrelaxed otherwise. Following this classification, they concluded that unrelaxed clusters are slightly worse modelled than relaxed ones, but that the mean bias of both types of clusters is completely consistent with that of the full sample.

5.6. Summary

In this section we have quantified the different contributions to the dispersion of mass estimates of clusters reconstructed from convergence maps. We find that M₅₀₀ cluster masses reconstructed from convergence maps have an uncertainty of the order of 9% that corresponds to the uncertainty in the fit and varies with the chosen model. In addition, a 10–16% uncertainty should be added to account for intrinsic projection effects. Regarding the variation of the bias along the cluster sample, which is of the order of 9% (${\sigma }_{b,\text{cluster-to-cluster}}^2$ = 0.009 for NFW and 0.007 for Einasto), we conclude that it is probably overestimated due to the limited projections per cluster available. When we try to correlate the bias and its dispersion to intrinsic characteristics of clusters, we find that redshift, mass, and dynamical state appear vaguely correlated to the bias. Additional effects that we have not considered in this work, such as the presence of substructures along the line of sight, could also be at the origin of the cluster-to-cluster scatter in observations.

The average bias of the mass estimates reconstructed with the method presented in this paper, as well as the fit uncertainty, are directly impacted by the radial range that is considered in the fit. Depending on the size of the error bars of the Σ-profile at each radius and the ability of the density model to describe the mass at the desired radii, mass estimates can be more or less biased. On the contrary, the intrinsic projection effect (quantified by ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$, as well as by the comparison of the biases obtained for the 0_pr_axes, 1_pr_axes, and 2_pr_axes projections) remains of the same order, even for the extreme case where we consider all the radial ranges available for the fit of the mass density model.

6. Comparison to gas observables

As investigated in several works (e.g. Meneghetti et al. 2010b; Rasia et al. 2012), the spatial distribution of the gas in the ICM does not follow necessarily the same distribution as the rest of the matter. This implies that the aforementioned conclusions may not apply to masses reconstructed from gas observables. In this section, we compare the projection effect obtained for total matter observables to the effect for gas observables.

To measure the impact of the projection when estimating masses from gas maps we take two approaches. The most direct option consists in measuring the mass of the gas from the gas mass maps (Sect. 3.2.3). We compute $M_{500}^{\mathrm{Gas},\mathrm{cyl}}$ by integrating the maps up to θ₅₀₀. The distribution of dispersions for the three masses calculated from the different projections is shown for all the clusters in the magenta histogram in Fig. 8. These results reflect the impact of the orientation as regards the arrangement of the gas in the cluster, but we note that M^Gas is not an estimate of its total mass.

A less direct option consists in estimating the mass from the integrated Y₅₀₀ signal in the y-maps by applying a Y₅₀₀ − M₅₀₀ scaling relation. First, the signal in the y-maps is integrated up to the real θ₅₀₀ for each cluster, which gives the cylindrically integrated Compton-y parameter $Y_{500}^{\mathrm{cyl}}$. We then convert the $Y_{500}^{\mathrm{cyl}}$ in M₅₀₀ assuming a given scaling relation.

There are several scaling relations in the literature that relate Y₅₀₀ to M₅₀₀ (A10, Cui et al. 2018). However, most of them look for a relation between the observed $Y_{500}^{\mathrm{cyl}}$ and the spherical mass $M_{500}^{\mathrm{sph}}$. In our case, we build a $Y_{500}^{\mathrm{cyl}}-M_{500}^{\kappa ,\mathrm{cyl}}$ scaling relation for our sample that allows us to compare cylindrically integrated masses. This scaling relation, together with $Y_{500}^{\mathrm{cyl}}-M_{500}^{\mathrm{sph}}$ are shown in Fig. E.1. The dispersion of the recovered masses for different projections of a cluster allows us to quantify the intrinsic projection effect. The empty green histogram in Fig. 8 shows the distribution of the dispersion of the biases estimated for all clusters from y-maps.

We observe that both approaches suggest a ${\sigma }_{b,\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}^2$ ∼ 10⁻⁴ (a few percent scatter on M₅₀₀) for the gas, meaning that for the considered data and within θ₅₀₀, the gas is an order of magnitude less disperse (i.e. more spherically distributed) than the dark matter also within the same θ₅₀₀. This is in agreement with previous works (Buote & Humphrey 2012; Becker & Kravtsov 2011; Pratt et al. 2019; Meneghetti et al. 2010b; Rasia et al. 2012) that also indicate that mass reconstructions from gas are less scattered than those from dark matter observables due to the difference in the spatial distribution of matter.

From the comparison between the bias of cylindrically integrated masses from the κ-maps and the bias of masses obtained from y-maps (with the $Y_{500}^{\mathrm{cyl}}-M_{500}^{\kappa ,\mathrm{cyl}}$ scaling relation), we obtain that their correlation is weak, with a Pearson correlation coefficient of 0.26. The correlation is stronger between the bias dispersions (i.e. σ_b) measured from the three random projections in κ- and y-maps: r = 0.48.

7. Summary and conclusions

We have studied the bias of galaxy cluster mass estimates reconstructed by fitting three-dimensional NFW and Einasto density models to projected mass density profiles obtained from convergence maps. We have performed the analysis making use of THE THREE HUNDREDGADGET-X hydrodynamical simulation clusters, selected to be representative of the NIKA2 LPSZ sample. All the results shown here were obtained with the 122 clusters of the three twin samples combined (Sect. 2.2). We checked that conclusions do not vary from twin sample to twin sample.

We decided to perform the fits of the mass density profiles accounting only for the radial ranges with uncertainties comparable to the ones of the profiles computed from CLASH convergence maps. Although the flexibility of the Einasto model permits fitting the projected density profiles by accounting for the different slopes, the fit fails describing the mass at ∼R₅₀₀, almost in the same way as the NFW model. Some works in the literature (Meneghetti et al. 2014) choose ad hoc the radial range and conclude that, as expected, models with more parameters fit better the density profile and therefore, give lower biases. The latter is also true for our case, but the difference between Einasto and NFW biases is not significant. Our main conclusion is that a projected mass density fit that is overall good does not give necessarily a good M₅₀₀ estimate.

Despite the slight differences between NFW and Einasto mass reconstructions, we have observed that both the mean bias and the scatter of the recovered masses are strongly correlated for the two models. This could mean that, in spite of the impact of the model choice, the effect is not enough to blur the information in the convergence maps. Quantitatively, the NFW and Einasto M₅₀₀ masses are biased by 8% and 7%, respectively. Therefore, both are in the 5 − 10% bias range obtained in previous simulation-based works (Rasia et al. 2012; Meneghetti et al. 2014; Becker & Kravtsov 2011; Euclid Collaboration 2024).

Regarding the errors associated with the mass estimates, we consider different contributions. With our approach, the uncertainty of the density profile fit to the convergence map introduces a mean error on the mass of 8% and 9% for NFW and Einasto profiles, respectively. An additional contribution comes from projection effects and we estimate that about 10–16% of dispersion should be considered. Such result matches the conclusions for lensing masses in the Pratt et al. (2019) review. Accounting for the projection effect together with the uncertainties of the fitting, the scatter is of the order of ∼18%, slightly less than the values in Becker & Kravtsov (2011) and Rasia et al. (2012), the latter being obtained from the fit of tangential shear profiles of simulated clusters.

When accounting for the full sample of THE THREE HUNDRED-NIKA2 LPSZ clusters, there is an excess of dispersion with respect to projection effects. We investigate the origin of such dispersion and we find, firstly, an uncertainty due to the fact that we use only three projections per cluster. Making use of 100 random projections for four clusters, we verify that considering only 3 projections per cluster, the mean bias of the three can be scattered about up to 7–10%.

We tried to correlate the bias to intrinsic properties of clusters and we checked the evolution of the mass bias with the true mass and redshift of clusters. Nevertheless, we did not find any evolutionary trend of the bias with the mass nor the redshift. In addition, we investigated the relation between the dynamical state and reconstructed masses of clusters, concluding that disturbed clusters are not particularly more biased than relaxed ones.

Regarding the orientation of clusters, we confirm that with our method and for clusters observed elongated along the line of sight, the reconstructed masses are overestimated. On the contrary, masses are underestimated if the major axes of clusters are on the plane of the sky. These conclusions are in agreement with the results in Euclid Collaboration (2024), where THE THREE HUNDRED convergence maps are used to create weak lensing observables for the preparation of the Euclid mission.

Given the similarities and differences between both works, it was indispensable to compare, all along this paper, our outcomes to Euclid Collaboration (2024). Despite all the presented nuances, both works seem to converge towards resembling conclusions.

Finally, we also compared how spherical is the spatial distribution of total matter and gas by measuring the dispersion from projection to projection in total mass and gas maps. Within R₅₀₀ the gas is more spherically distributed than the dark matter, which allows one to have a mass reconstruction from projected maps less dependent on the orientation of the cluster. This was already known from previous studies (Meneghetti et al. 2010b), which motivated the observations of clusters in X-rays and SZ. For observational HSE-to-lensing mass biases estimated as for the LPSZ clusters (Sect. 2.1), as well as for other works in the literature (e.g. Bartalucci et al. 2018; Sereno & Ettori 2015; Foëx et al. 2012), this would imply that most of the projection scatter is introduced by the lensing mass estimate.

Acknowledgments

This work has been made possible by THE THREE HUNDRED collaboration (https://the300-project.org). M.D.P., A.F., and A.P. acknowledge support from Sapienza Università di Roma thanks to Progetti di Ricerca Medi 2021, RM12117A51D5269B. W.C. is supported by the STFC AGP Grant ST/V000594/1 and the Atracción de Talento Contract no. 2020-T1/TIC-19882 granted by the Comunidad de Madrid in Spain. W.C. and G.Y. thank the Ministerio de Ciencia e Innovación (Spain) for financial support under Project grant PID2021-122603NB-C21. This project was carried out using the Python libraries matplotlib (Hunter 2007), numpy (Harris et al. 2020), astropy (Astropy Collaboration 2013, 2018), and emcee (Foreman-Mackey et al. 2019; Goodman & Weare 2010).

References

Appendix A: Mass reconstruction

A.1. Mass density model best-fit parameters

In Fig. A.1 we present the relation between the best-fit values of the posteriors for the parameters of all NFW and Einasto fits and we give the Pearson correlation coefficient for each pair of parameters. We define the best-fit parameters from the maximum point of multi-dimensional posterior distributions. In the left panel, we show with big stars the mean concentration per M₂₀₀ bin for NFW and compare to the concentration-mass relations in Cui et al. (2018). The concentration-mass relation obtained in this work with NFW follows the flat relation of the GADGET-X clusters results from THE THREE HUNDRED in Cui et al. (2018). See Sect. 4.1 for details. As a reference, we also show the horizontal c₂₀₀ = 3 line that corresponds to the concentration value fixed in Euclid Collaboration (2024) to obtain less biased masses when fitting tNFW density models.

A.2. Reconstructed mass profiles

In Fig. A.2 we show the mass profiles corresponding to the projected mass density fits of Fig. 3 with NFW and Einasto models for the 0306 cluster in the snapshot 101. We present the profiles for the projections along the 3 main axes of the simulation on top, and in the bottom the profiles for the main inertia moment axes of the cluster. The dashed line shows the R₅₀₀ − M₅₀₀ relation given by the definition: $M_{500}/\frac{4}{3}\pi R_{500}^3=500{\rho }_{\mathrm{crit}}$. The cyan profile shows the spherical true mass profile (Sect. 3.1). The departure of the NFW and Einasto profiles from the real spherical mass profile at ∼R₅₀₀ will determine the difference between the true M₅₀₀ and the NFW and Einasto estimates.

Fig. A.1. Best-fit values of the posterior distributions of NFW (left) and Einasto (right) parameters. Each point corresponds to one projection of one cluster. We give the Pearson correlation coefficients for the parameters presented in each panel. In the left panel the big stars show the mean concentration and the standard deviation in different mass bins. Blue and red lines show the concentration-mass relations for GADGET-X and GADGET-MUSIC simulations obtained in Cui et al. (2018) accounting for different radial ranges: with solid lines considering data above 0.05 × R200 and with dashed lines the results considering data above 34h−1 kpc. The horizontal orange line at c200 = 3 in the left panel represents the best concentration value obtained in Euclid Collaboration (2024) to reconstruct unbiased masses. For comparison, in the left panel we also plot the M200 and c200 = R200/rs obtained for the Einasto best-fits.

Fig. A.2. Mass profiles reconstructed from the convergence map fits for the 0306 cluster in the snapshot 101 (z = 0.817). The orange and purple profiles correspond to NFW and Einasto models, respectively. We give the mean profiles with 1σ contours. The cyan profile is the spherical mass profile obtained from the simulation. The black dashed line shows the R500 − M500 relation.

Appendix

Mass bias scatter distributions

We present in Fig. B.1 the distributions of ${\sigma }_{b,\mathrm{f}\mathrm{i}\mathrm{t}}^2$ (left) and ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ (right) for all the clusters in our sample. The histograms correspond to the bias variance for the NFW (orange) and Einasto (purple) models.

Fig. B.1. Variance of mass bias for all clusters. Left: scatter induced by uncertainties on the fit, for all the projections of all the clusters in our sample. Right: variance of the bias in the mass reconstruction of each cluster.

Appendix C: Validation of the projection effect

As we are interested in checking the real scatter due to the projection effects, we produced for 4 different clusters 100 κ-maps along random projections. We chose one cluster per snapshot, with diverse dynamical states, mostly unrelaxed (Table C.1, dynamical state indicators defined following the definitions in De Luca et al. 2021).

For each map we followed the NFW and Einasto mass reconstruction procedure described in Sect. 4.1. Therefore, for the considered four clusters, we measured the κ-bias of each projection. The mean bias and variance for the 100 projections, for the 3 main axes and for the 3 inertia moment axes are summarised in Table C.2. Although considering the principal axes of the clusters regarding their moments of inertia is a way to account for the whole dispersion due to projection, it gives overestimated scatter values with respect to the dispersion obtained from 100 projections. The dispersion, when estimated with only 3 data points, can be very unreliable, but for the 3 inertia moment axes it tends to be overestimated, while for the 3 main (that is, random) axes it can be both over and underestimated. Hence, for the analysis with the full sample of clusters we decided to remove the 0_pr_axes, 1_pr_axes, and 2_pr_axes projections.

Appendix D: Scatter from spherical mass estimates and from projected masses within a fixed aperture

In this section we demonstrate that masses cylindrically integrated within a fixed aperture are less scattered than spherically integrated masses. We use simulated mass density profiles of different shapes following the NFW (Eq. 5), gNFW, tNFW, and Hernquist (Hernquist 1990) models.

We define the gNFW profile as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{gNFW}}(r)=\frac{{\rho }_s}{{\left(\frac{r}{r_s}\right)}^{\gamma }{(1+{\left(\frac{r}{r_s}\right)}^{\alpha })}^{\frac{\beta -\gamma }{\alpha }}},\end{array}$$(D.1)

with α = 3, r_s = R_Δ/c_Δ and:

$$\begin{array}{c}\hfill {\rho }_s=\frac{c_{\mathrm{\Delta }}^3\mathrm{\Delta }{\rho }_{\mathrm{crit}}(-3+\gamma )}{-3{(c_{\mathrm{\Delta }}^3)}^{1-\gamma /3}{}_2{F}_1((\beta -\gamma )/3,1-\gamma /3,2-\gamma /3,-c_{\mathrm{\Delta }}^3)}.\end{array}$$(D.2)

Here ₂F₁ is a Gaussian hypergeometric function. The tNFW is given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{tNFW}}(r)=\frac{{\rho }_s}{r/r_s{(1+r/r_s)}^2}{\left(\frac{{\tau }^2}{{\tau }^2+{(r/r_s)}^2}\right)}^{\eta },\end{array}$$(D.3)

with η = 1, r_s = R_Δ/c_Δ and:

$$\begin{array}{c}\hfill {\rho }_s=\frac{c_{\mathrm{\Delta }}^3\mathrm{\Delta }{\rho }_{\mathrm{crit}}}{3I(c_{\mathrm{\Delta }})},\end{array}$$(D.4)

where,

$$\begin{array}{cc}\hfill I(c_{\mathrm{\Delta }})& ={\tau }^2\frac{-1[2(1+{\tau }^2)-(-1+{\tau }^2)ln({\tau }^2)]}{2{(1+{\tau }^2)}^2}+\hfill \\ \hfill & {\tau }^2\frac{\frac{2(1+{\tau }^2)}{(1+c_{\mathrm{\Delta }})}+4\tau \arctan (c_{\mathrm{\Delta }}/\tau )}{2{(1+{\tau }^2)}^2}+\hfill \\ \hfill & {\tau }^2\frac{(-1+{\tau }^2)[2ln(1+c_{\mathrm{\Delta }})-ln(c_{\mathrm{\Delta }}^2+{\tau }^2)]}{2{(1+{\tau }^2)}^2}.\hfill \end{array}$$(D.5)

The Hernquist mass density profile is defined as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{Hernquist}}(r)=\frac{{\rho }_s}{\frac{r}{r_s}{(1+\frac{r}{r_s})}^3},\end{array}$$(D.6)

with r_s = R_Δ/c_Δ and:

$$\begin{array}{c}\hfill {\rho }_s=\frac{c_{\mathrm{\Delta }}\mathrm{\Delta }{\rho }_{\mathrm{crit}}2{(1+c_{\mathrm{\Delta }})}^2}{3}.\end{array}$$(D.7)

We produce profiles with a wide variety of parameters for 5 different initial M₅₀₀ values between 3 × 10¹⁴ M_⊙ and 11 × 10¹⁴ M_⊙ and for the four different redshifts in our twin samples (z = [0.490, 0.592, 0.700, 0.817]). For NFW and Hernquist models we vary c₅₀₀ from 1 to 6. For gNFW and tNFW we fix c₅₀₀ = 4 and vary β = 3 − 6 and γ = 0.4 − 1 for gNFW and τ = 0.2 − 4 for tNFW. All the considered projected mass density profiles are shown in the left panel in Fig. D.1 to demonstrate the variety of the analysed density shapes.

To mimic a departure from the true mass, for each profile with M₅₀₀ we estimate a distribution of new profiles for masses that vary from −18% to 18% from the true M₅₀₀, $M_{500}^{\mathrm{new}\mathrm{sph}}$ (based on an average ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ that we measured from data). For each of the new profiles we can integrate the projected density up to the original θ₅₀₀, $M_{500}^{\mathrm{new}\mathrm{cyl}}$. Comparing the new $M_{500}^{\mathrm{new}\mathrm{sph}}$ and $M_{500}^{\mathrm{new}\mathrm{cyl}}$ to the original M₅₀₀ we get the spherical and cylindrical biases. From the distribution of biased masses, we take randomly trios to simulate the 3 random projections and estimate their variance: ${\sigma }_{b^{\mathrm{s}\mathrm{p}\mathrm{h}}}^2$ is the variance of three biased spherical masses and with their corresponding cylindrically integrated masses we get ${\sigma }_{b^{\mathrm{c}\mathrm{y}\mathrm{l}}}^2$.

In the right panel in Fig. D.1 we present the relation between the mean ${\sigma }_{b^{\mathrm{s}\mathrm{p}\mathrm{h}}}^2$ and ${\sigma }_{b^{\mathrm{c}\mathrm{y}\mathrm{l}}}^2$ for each mass, redshift, and profile shape considered. Same colours as in Fig. D.1 are used to make reference to each density profile type. As expected, the figure shows that the dispersion of cylindrically integrated masses at a given aperture tends to be smaller (${\sigma }_{b^{\mathrm{s}\mathrm{p}\mathrm{h}}}^2-{\sigma }_{b^{\mathrm{c}\mathrm{y}\mathrm{l}}}^2\sim 0.001$ to 0.003) than spherical integrations up to the corresponding radius in each case. The blue points close to the one-to-one line correspond to NFW models with a concentration of c₅₀₀ = 1.

Fig. D.1. Projected mass density profiles used to quantify the difference in scatter for spherically and cylindrically integrated masses (left). The variance of cylindrically integrated masses with respect to spherically integrated masses (right). The black dashed line corresponds to ${\sigma }_{b^{\mathrm{s}\mathrm{p}\mathrm{h}}}^2={\sigma }_{b^{\mathrm{c}\mathrm{y}\mathrm{l}}}^2$. In blue NFW, in red gNFW, in green tNFW, and in orange Hernquist models.

Appendix E: Building the Y₅₀₀ − M₅₀₀ scaling relation

We construct our own scaling relation following the parametrisation in Eq. 16 in Arnaud et al. (2010):

$$\begin{array}{c}\hfill h{(z)}^{-2/3}{Y}_{500}=\beta {\left[\frac{M_{500}}{3\times {10}^{14}{h}_{70}^{-1}{\mathrm{M}}_{\odot }}\right]}^{\alpha }{h}_{70}^{-1},\end{array}$$(E.1)

where Y₅₀₀ is in Mpc² and M₅₀₀ in M_⊙. Here h(z) is the ratio of the Hubble constant at redshift z to its present value, and, h₇₀ the Hubble constant in units of 70 km/s/Mpc.

We compare two scaling relations that relate the cylindrical or projected Y₅₀₀ to the spherical and cylindrical M₅₀₀ mass, respectively. The spherical quantities are computed from the ICM profiles: $M_{500}^{\mathrm{sph}}$ is the mass profile evaluated at R₅₀₀, that is, the true M₅₀₀ throughout this paper. The cylindrical quantities $M_{500}^{\mathrm{cyl}}$ and $Y_{500}^{\mathrm{cyl}}$ are obtained by integrating the signal within θ₅₀₀ in the κ- and y-maps.

In each case, we perform a fit of Eq. E.1 using the orthogonal BCES method (Akritas & Bershady 1996). Uncertainties on fit parameters are estimated from 100000 bootstrap resamples of the data. We show the scaling relations in Fig. E.1, as well as the best-fit parameters (top legend in the figure). The slopes of both relations are consistent, but since the cylindrically integrated masses are larger than spherically integrated ones, the intercepts are not compatible.

We compare, as in Fig. 10 in Cui et al. (2018), the $Y_{500}^{\mathrm{cyl}}-M_{500}^{\mathrm{sph}}$ relation to results in the literature. In Cui et al. (2018), the authors computed the scaling relation using all THE THREE HUNDRED clusters at z = 0, for both GADGET-X and GADGET-MUSIC flavours. The scaling relation obtained with our sample has a lower normalisation and slightly steeper slope than the results obtained with GADGET-X simulations in Cui et al. (2018). The differences may be due to both the considered cluster sample and the redshift range. The GADGET-MUSIC scaling relation from Cui et al. (2018) has an even smaller intercept value, showing the impact of the hydrodynamical model on the simulation. In the same figure we also compare the results to the scaling relations obtained from observations in Planck Collaboration XX (2014) and Nagarajan et al. (2018). Our scaling relation lies between both observational results, with a flatter slope than Planck Collaboration XX (2014) and steeper than Nagarajan et al. (2018). Given the offset between spherical $M_{500}^{\mathrm{sph}}$ and cylindrical $M_{500}^{\mathrm{cyl}}$ masses, the observational results shown in the right panel in Fig. E.1 are shifted with respect to the scaling relation in the left panel and are not comparable.

Fig. E.1. Scaling relations between Y500 and M500 for cylindrically integrated SZ signal with respect to spherically (right) and cylindrically (left) integrated masses. Black dots correspond to the values per cluster in our sample, while the blue line shows the best-fit relation. The solid brown line shows the result in Planck Collaboration XX (2014), the black dash-dotted is from Nagarajan et al. (2018). The pink dotted and green dashed relations are the results in Cui et al. (2018) for the GADGET-MUSIC and GADGET-X simulations.

All Tables

All Figures

Fig. 1. Normalised gas pressure (left) and total matter mass density (right) spherical profiles for all the clusters in THE THREE HUNDRED-NIKA2 LPSZ twin samples. We also plot the “universal” pressure profile from the X-ray analysis in A10 and the X-ray and SZ analysis in P13. All profiles are drawn with respect to the radii normalised by the R500 corresponding to each cluster. The grey shaded areas indicate the regions below 0.025 × R500 that are not numerically robust.

In the text

Fig. 2. Variance of the convergence profiles (left) and fractional error of the Σ-profiles (right) calculated from convergence maps. In grey we show the profiles for the simulated clusters in THE THREE HUNDRED-NIKA2 LPSZ samples. Red, green, blue, and yellow profiles correspond to CL J1226.9+3332, PSZ2 G144.83+25.11, PSZ2 G228.16+75.20, and MACS J1423.8+2404 clusters, respectively, with data from the CLASH convergence maps (Zitrin et al. 2015) and uncertainties computed following Muñoz-Echeverría et al. (2023; sample defined in Muñoz-Echeverría et al. 2022). Solid and dashed lines differentiate the profiles for the CLASH convergence maps reconstructed assuming the Light-Traces-Mass (LTM) and the Pseudo-Isothermal Elliptical Mass Distribution with an elliptical NFW (PIEMD+eNFW) models, respectively. Vertical dotted lines in the right panel indicate the minimum and maximum of the used radial range.

In the text

Fig. 3. Projected mass density profiles (empty markers) for the six convergence maps for the 0306 cluster in the snapshot 101 (z = 0.817). The orange and purple profiles show the fit with NFW and Einasto models, respectively. We give the mean profiles with 1σ contours. From left to right and top to bottom we present the 0, 1, and 2 random projections and 0_pr_axes, 1_pr_axes, and 2_pr_axes principal moment of intertia projections. The vertical cyan line shows the true R500 of the cluster.

In the text

Fig. 4. Mass bias for masses inferred from fitted κ-maps with respect to M500 per cluster. Each data point corresponds to the mean bias of a cluster for a given random projection (0, 1 or 2) with error bars showing the 1σ uncertainty. Results for the Einasto model have been artificially shifted in mass for visualisation purposes. The horizontal orange and purple dashed lines show the mean for all the results for NFW and Einasto, respectively. The black line indicates the zero.

In the text

	Fig. 5. Relative differences of Σ-profiles between the fits with NFW (top) and Einasto (bottom) models and the profiles measured from the convergence maps. Solid lines show the median.
In the text

	Fig. 6. Mean bias of the three projections (left) and ${\sigma }_{b,\text{3}\hspace{0.17em}\text{projections}}^2$ (right) per cluster for Einasto against NFW reconstructions. Each dot corresponds to one cluster. We give the Pearson correlation coefficient for the values given in each panel.
In the text

Fig. 7. Mass bias for NFW (orange) and Einasto (purple) fits of the Σ-profiles when aligned with the principal axes of the clusters regarding their moments of inertia. Solid lines correspond to 0_pr_axes, dashed lines to 1_pr_axes, and dotted lines to 2_pr_axes. We give the value of the median and the 16th and 84th percentiles. Shaded areas represent the 1σ regions of the bias obtained for all the clusters with the three random projections (Fig. 4).

In the text

Fig. 8. Variance of the mass bias estimated from maps. The grey histogram is obtained from the integration of convergence maps, i.e. $M_{500}^{\kappa ,\mathrm{cyl}}$. The magenta histogram corresponds to the variance for gas masses ($M_{500}^{\mathrm{Gas},\mathrm{cyl}}$) estimated from gas mass maps. The empty green histogram shows the dispersion of the $M_{500}^{\mathrm{cyl}}$ after conversion from the $Y_{500}^{\mathrm{cyl}}$ measured from y-maps, using the scaling relation in Fig. E.1.

In the text

	Fig. 9. Mean bias per cluster. The vertical lines show the mean bias when using 100 projections to reconstruct the mass. The histograms are obtained from the mean bias computed by taking randomly 3 projections out of the 100 available 100 times.
In the text

	Fig. 10. Bias for clusters in different redshift (left) and mass (right) bins. We give the mean value and the σb, mean scatter for the mean bias of clusters in each bin. Results for NFW and Einasto models are shown in orange and purple, respectively, with a line as a guide to the eye.
In the text

Fig. 11. Bias properties depending on the dynamical state of clusters within R200. Top: distributions of the mean bias per cluster for relaxed (green), intermediate (blue), and disturbed (red) clusters. We show the bias from NFW (left) and Einasto (right) reconstructions. Bottom: distributions of the variance of the bias per cluster accounting for three projections. Clusters are classified according to their dynamical state and results for NFW (left) and Einasto (right) are shown. The vertical lines show the mean values of the distributions.

In the text

Fig. 12. Summarised bias properties depending on the dynamical state of clusters within R200 (empty) and R500 (filled). Left: mean (circle) and median (star) of the mean biases of each type of cluster. Right: mean and median of the ${\sigma }_{b,3\mathrm{projections}}^2$ for all the clusters of each type. We show in green, blue, and red the results for relaxed, intermediate, and disturbed clusters, respectively.

In the text

Fig. A.1. Best-fit values of the posterior distributions of NFW (left) and Einasto (right) parameters. Each point corresponds to one projection of one cluster. We give the Pearson correlation coefficients for the parameters presented in each panel. In the left panel the big stars show the mean concentration and the standard deviation in different mass bins. Blue and red lines show the concentration-mass relations for GADGET-X and GADGET-MUSIC simulations obtained in Cui et al. (2018) accounting for different radial ranges: with solid lines considering data above 0.05 × R200 and with dashed lines the results considering data above 34h−1 kpc. The horizontal orange line at c200 = 3 in the left panel represents the best concentration value obtained in Euclid Collaboration (2024) to reconstruct unbiased masses. For comparison, in the left panel we also plot the M200 and c200 = R200/rs obtained for the Einasto best-fits.

In the text

	Fig. A.2. Mass profiles reconstructed from the convergence map fits for the 0306 cluster in the snapshot 101 (z = 0.817). The orange and purple profiles correspond to NFW and Einasto models, respectively. We give the mean profiles with 1σ contours. The cyan profile is the spherical mass profile obtained from the simulation. The black dashed line shows the R500 − M500 relation.
In the text

	Fig. B.1. Variance of mass bias for all clusters. Left: scatter induced by uncertainties on the fit, for all the projections of all the clusters in our sample. Right: variance of the bias in the mass reconstruction of each cluster.
In the text

Fig. D.1. Projected mass density profiles used to quantify the difference in scatter for spherically and cylindrically integrated masses (left). The variance of cylindrically integrated masses with respect to spherically integrated masses (right). The black dashed line corresponds to ${\sigma }_{b^{\mathrm{s}\mathrm{p}\mathrm{h}}}^2={\sigma }_{b^{\mathrm{c}\mathrm{y}\mathrm{l}}}^2$. In blue NFW, in red gNFW, in green tNFW, and in orange Hernquist models.

In the text

Fig. E.1. Scaling relations between Y500 and M500 for cylindrically integrated SZ signal with respect to spherically (right) and cylindrically (left) integrated masses. Black dots correspond to the values per cluster in our sample, while the blue line shows the best-fit relation. The solid brown line shows the result in Planck Collaboration XX (2014), the black dash-dotted is from Nagarajan et al. (2018). The pink dotted and green dashed relations are the results in Cui et al. (2018) for the GADGET-MUSIC and GADGET-X simulations.

In the text