# Exploring the hydrostatic mass bias in MUSIC clusters: application to the NIKA2 mock sample

Giulia Gianfagna, Marco de Petris, Gustavo Yepes, Federico de Luca, Federico

Sembolini, Weiguang Cui, Veronica Biffi, Florian Kéruzoré, Juan

Macías-Pérez, Frédéric Mayet, et al.

► To cite this version:

Giulia Gianfagna, Marco de Petris, Gustavo Yepes, Federico de Luca, Federico Sembolini, et al. Exploring the hydrostatic mass bias in MUSIC clusters: application to the NIKA2 mock sample. Monthly Notices of the Royal Astronomical Society, 2021, 502 (4), pp.5115-5133. 10.1093/mnras/stab308. hal-02973179

HAL Id: hal-02973179

https://hal.science/hal-02973179

Submitted on 20 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Exploring the hydrostatic mass bias in MUSIC clusters: application to the NIKA2 mock sample

Giulia Gianfagna,<fnref n="1" />* Marco De Petris,<fnref n="1" /> Gustavo Yepes,<fnref n="2" /> Federico De Luca,<fnref n="1" />,<fnref n="3" /> Federico Sembolini,<fnref n="1" />,<fnref n="2" />

Weiguang Cui,<fnref n="4" /> Veronica Biffi,<fnref n="5" />,<fnref n="6" /> Florian Kéruzoré,<fnref n="7" /> Juan Macías-Pérez,<fnref n="7" /> Frédéric Mayet,<fnref n="7" />
Laurence Perotto,<fnref n="7" /> Elena Rasia<fnref n="6" />,<fnref n="8" /> and Florian Ruppin<fnref n="9" />

Accepted 2021 February 1. Received 2021 January 18; in original form 2020 October 7

## ABSTRACT

Clusters of galaxies are useful tools to constrain cosmological parameters, only if their masses can be correctly inferred from observations. In particular, X-ray and Sunyaev–Zeldovich (SZ) effect observations can be used to derive masses within the framework of the hydrostatic equilibrium. Therefore, it is crucial to have a good control of the possible mass biases that can be introduced when this hypothesis is not valid. In this work, we analysed a set of 260 synthetic clusters from the MUSIC simulation project at redshifts 0 ≤ *z* ≤ 0.82. We estimate the hydrostatic mass of the MUSIC clusters from X-ray only (temperature and density) and from X-ray and SZ (density and pressure). Then, we compare them with the true 3D dynamical mass. The biases are of the order of 20 per cent. We find that using the temperature instead of the pressure leads to a smaller bias, although the two values are compatible within 1σ. Non-thermal contributions to the total pressure support, arising from bulk motion and turbulence of the gas, are also computed and show that they are sufficient to account for this bias. We also present a study of the correlation between the mass bias and the dynamical state of the clusters. A clear correlation is shown between the relaxation state of the clusters and the bias factor. We applied the same analysis on a subsample of 32 objects, already selected for supporting the NIKA2 SZ Large Program.

**Key words:** methods: numerical – galaxies: clusters: general – galaxies: clusters: intracluster medium – large-scale structure of Universe.

## 1 INTRODUCTION

Galaxy clusters are the most massive gravitationally bound objects in the Universe and they are mainly composed by dark matter, that amounts to 80 per cent of the total mass (for a full review, see Kravtsov & Borgani 2012). About 8 per cent is composed by galaxies and the remaining 12 per cent is represented by the so-called intracluster medium (ICM), i.e. the hot gas located between galaxies. This gas component provides significant physical information, as it can be observed in the X-ray band and in millimetre wavelengths through the Sunyaev–Zeldovich (SZ) effect (Sunyaev & Zeldovich 1972).

The emission in the X-ray band is mainly due to the thermal bremsstrahlung. From this emission, we can directly measure the temperature, which determines the bremsstrahlung cut-off, and the

electron density, to which the spectrum normalization is proportional (for a review, see Boehringer & Werner 2009). X-ray observations occurred to be particularly successful because the emission is proportional to the square of the gas density (for a review on the methods adopted to reconstruct the mass profiles in X-ray luminous galaxy clusters, see Ettori et al. 2013). In the last two decades, X-ray observatories with improved sensitivity and angular resolution, like *XMM-Newton* and *Chandra*, has conducted cluster X-ray emission studies over large areas of the sky (Pacaud et al. 2007; Andreon et al. 2017), and deeper studies of previously known objects (e.g. Vikhlinin et al. 2009; Mantz et al. 2010).

X-ray observations are mainly exploring the central regions of the clusters, where the electron number density is high, although there are also X-ray projects exploring the cluster outskirts, such as X-COP (Eckert et al. 2017). A more efficient way to map the cluster outskirts is through the thermal SZ effect (Sunyaev & Zeldovich 1972). This effect depends linearly on the pressure, so it is less sensitive to the density decrease at high radial distances from the cluster centre. The

<footnote n="1">¹Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro, I-5-00185 Roma, Italy</footnote>

<footnote n="2">²Departamento de Física Teórica and CIAFF, Módulo 8, Facultad de Ciencias, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="3">³Dipartimento di Fisica, Università di Roma 'Tor Vergata', Via della Ricerca Scientifica, I-00133 Roma, Italy</footnote>

<footnote n="4">⁴Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="5">⁵Fakultät für Physik, Universität-Sternwarte München, LMU Munich, Scheinerstr 1, D-81679 München, Germany</footnote>

<footnote n="6">⁶IFPU – Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

<footnote n="7">⁷Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, avenue des Martyrs, F-38000 Grenoble, France</footnote>

<footnote n="8">⁸INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="9">⁹Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote>* E-mail: gianfagna.1665033@student.uniroma1.it</footnote>
<figure><img src="image_2.png" /></figure>

# Exploring the hydrostatic mass bias in MUSIC clusters: application to the NIKA2 mock sample

Giulia Gianfagna,<fnref n="1" />* Marco De Petris,<fnref n="1" /> Gustavo Yepes,<fnref n="2" /> Federico De Luca,<fnref n="1" />,<fnref n="3" /> Federico Sembolini,<fnref n="1" />,<fnref n="2" />
Weiguang Cui,<fnref n="4" /> Veronica Biffi,<fnref n="5" />,<fnref n="6" /> Florian Kéruzoré,<fnref n="7" /> Juan Macías-Pérez,<fnref n="7" /> Frédéric Mayet,<fnref n="7" />
Laurence Perotto,<fnref n="7" /> Elena Rasia<fnref n="6" />,<fnref n="8" /> and Florian Ruppin<fnref n="9" />

Accepted 2021 February 1. Received 2021 January 18; in original form 2020 October 7

## ABSTRACT

Clusters of galaxies are useful tools to constrain cosmological parameters, only if their masses can be correctly inferred from observations. In particular, X-ray and Sunyaev–Zeldovich (SZ) effect observations can be used to derive masses within the framework of the hydrostatic equilibrium. Therefore, it is crucial to have a good control of the possible mass biases that can be introduced when this hypothesis is not valid. In this work, we analysed a set of 260 synthetic clusters from the MUSIC simulation project at redshifts 0 ≤ z ≤ 0.82. We estimate the hydrostatic mass of the MUSIC clusters from X-ray only (temperature and density) and from X-ray and SZ (density and pressure). Then, we compare them with the true 3D dynamical mass. The biases are of the order of 20 per cent. We find that using the temperature instead of the pressure leads to a smaller bias, although the two values are compatible within 1σ. Non-thermal contributions to the total pressure support, arising from bulk motion and turbulence of the gas, are also computed and show that they are sufficient to account for this bias. We also present a study of the correlation between the mass bias and the dynamical state of the clusters. A clear correlation is shown between the relaxation state of the clusters and the bias factor. We applied the same analysis on a subsample of 32 objects, already selected for supporting the NIKA2 SZ Large Program.

**Key words:** methods: numerical – galaxies: clusters: general – galaxies: clusters: intracluster medium – large-scale structure of Universe.

## 1 INTRODUCTION

Galaxy clusters are the most massive gravitationally bound objects in the Universe and they are mainly composed by dark matter, that amounts to 80 per cent of the total mass (for a full review, see Kravtsov & Borgani 2012). About 8 per cent is composed by galaxies and the remaining 12 per cent is represented by the so-called intracluster medium (ICM), i.e. the hot gas located between galaxies. This gas component provides significant physical information, as it can be observed in the X-ray band and in millimetre wavelengths through the Sunyaev–Zeldovich (SZ) effect (Sunyaev & Zeldovich 1972).

The emission in the X-ray band is mainly due to the thermal bremsstrahlung. From this emission, we can directly measure the temperature, which determines the bremsstrahlung cut-off, and the

electron density, to which the spectrum normalization is proportional (for a review, see Boehringer & Werner 2009). X-ray observations occurred to be particularly successful because the emission is proportional to the square of the gas density (for a review on the methods adopted to reconstruct the mass profiles in X-ray luminous galaxy clusters, see Ettori et al. 2013). In the last two decades, X-ray observatories with improved sensitivity and angular resolution, like *XMM-Newton* and *Chandra*, has conducted cluster X-ray emission studies over large areas of the sky (Pacaud et al. 2007; Andreon et al. 2017), and deeper studies of previously known objects (e.g. Vikhlinin et al. 2009; Mantz et al. 2010).

X-ray observations are mainly exploring the central regions of the clusters, where the electron number density is high, although there are also X-ray projects exploring the cluster outskirts, such as X-COP (Eckert et al. 2017). A more efficient way to map the cluster outskirts is through the thermal SZ effect (Sunyaev & Zeldovich 1972). This effect depends linearly on the pressure, so it is less sensitive to the density decrease at high radial distances from the cluster centre. The

<footnote n="1">¹Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro, I-5-00185 Roma, Italy</footnote>

<footnote n="2">²Departamento de Física Teórica and CIAFF, Módulo 8, Facultad de Ciencias, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="3">³Dipartimento di Fisica, Università di Roma 'Tor Vergata', Via della Ricerca Scientifica, I-00133 Roma, Italy</footnote>

<footnote n="4">⁴Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="5">⁵Fakultät für Physik, Universität-Sternwarte München, LMU Munich, Scheinerstr 1, D-81679 München, Germany</footnote>

<footnote n="6">⁶IFPU – Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

<footnote n="7">⁷Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, avenue des Martyrs, F-38000 Grenoble, France</footnote>

<footnote n="8">⁸INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="9">⁹Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote>* E-mail: gianfagna.1665033@student.uniroma1.it</footnote>
SZ effect is a redshift independent probe that consists of the spectral distortion of the cosmic microwave background (CMB) radiation due to the inverse Compton scattering of the CMB photons with the hot electrons of the ICM (for a review, see Carlstrom, Holder & Reese 2002; Kitayama 2014; Mroczkowski et al. 2019). It can be used to directly measure the ICM pressure distribution. The latter can be combined with the electron density from an X-ray observation to infer the cluster mass profile. This method avoids the necessity of obtaining the deep X-ray observations required to measure a spatially resolved temperature profile. To date, the SZ effect induced by thermal electrons (tSZ) has been detected for more than a thousand galaxy clusters, including more than 200 new clusters previously unknown by any other observational means (Kitayama 2014) thanks to new observing facilities such as the South Pole Telescope (SPT; Williamson et al. 2011; Reichardt et al. 2013; Bleem et al. 2020), the Atacama Cosmology Telescope (ACT; Marriage et al. 2010; Hasselfield et al. 2013; Hilton et al. 2018, 2020), and the Planck satellite (Planck Collaboration VIII, XXIX, XXVII 2011a, 2014c, 2016b).

Both X-ray and SZ observations can be used to infer the mass of a cluster. In particular, from the former we can exploit the temperature and the electron density, from the latter the pressure profile. In order to use this information, we need to make three fundamental assumptions: the gas must trace the cluster potential well, it must be spherically symmetric, and in hydrostatic equilibrium (HE; Kravtsov & Borgani 2012 for a review). The mass inferred through this method is called hydrostatic mass (see Section 2). An easy way to track the error made when using HE is the hydrostatic mass bias *b* (Planck Collaboration XX 2014b). It is defined as the difference of the cluster total mass to the one estimated by HE, divided by the total mass *b* = ($M_{tot} - M_{HE}$)/$M_{tot}$.

Galaxy clusters, and, in particular, their number counts are a fundamental cosmological tool. The abundance of clusters and its evolution with redshift are particularly sensitive to the cosmic matter density, Ω<sub>m</sub>, and the present amplitude of density fluctuations, characterized by σ<sub>8</sub>, i.e. the rms linear overdensity in spheres of radius 8 h<sup>-1</sup> Mpc. The CMB primary anisotropies, on the other hand, are related to the density perturbation power spectrum at the time of recombination. A comparison of the amplitude of density perturbations from recombination until today, allows us to look for possible extensions to the concordance ΛCDM model, such as non-minimal neutrino masses or non-zero curvature contributions (Planck Collaboration XXIV 2016a; Salvati, Douspis & Aghanim 2018).

In Planck Collaboration XX (2014b), a tension between the number of clusters detected by SZ signal, and the number of clusters predicted from the cosmological parameters inferred from the primary CMB power spectrum is reported. This issue was later confirmed in Planck Collaboration XXIV (2016a) and arises when the HE mass bias is fixed to a constant value of 0.2 in cluster cosmological analyses. The hydrostatic mass bias *b* plays an important role in the number counts, because it leads to a modification of the cluster population at a given mass. In particular, it significantly affects the value of σ<sub>8</sub>: the lower the (1 − *b*), the higher the σ<sub>8</sub>, (Salvati et al. 2018; Ruppin et al. 2019a). Salvati et al. (2018) published an update of the constraints on cosmological parameters from the clusters observed by *Planck*. They find that the bias needed to reconcile CMB constraints with those from the tSZ number counts is (1 − *b*) = 0.62 ± 0.07, which is compatible with the value (1 − *b*) = 0.58 ± 0.04 found in Planck Collaboration XXIV (2016a). This value is confirmed by Koukoufilippas et al. (2020), who cross-correlate *Planck* maps of the tSZ Compton-y parameter with the galaxy distribution.

These values derived from observations are nevertheless in disagreement with the value of (1 − *b*) estimated from simulations, which is about 0.8. This topic will be deeply discussed later in this work. Moreover, there is a factor of 2.5 more clusters predicted than observed when taking into account the CMB cosmology and a value (1 − *b*) of 0.8 (Andreon 2014; Planck Collaboration XXIV 2016a; Salvati et al. 2018). This value of the bias, derived from simulations, has been recently confirmed by Makiya, Hikage & Komatsu (2020). They perform a joint analysis of power spectra of the tSZ and the cosmic weak lensing shear, in order to obtain a (1 − *b*) constraint which is independent from the primordial CMB spectrum. They find (1 − *b*) = 0.73<sup>+0.08</sup><sub>-0.13</sub>, and conclude that the late-time probe (tSZ and cosmic weak lensing shear) cosmologies are consistent with each other, but they could not be totally consistent with the CMB cosmology, which is leading to a different value of the mass bias.

However, the cluster number counts are limited by systematic effects, in particular those affecting the mass estimates. The tSZ power spectrum, in turn, is not measured with sufficient accuracy, especially at small angular scales, to reduce the tension with the CMB. The tSZ cosmological analysis can be improved by considering more realistic and complex hypotheses on the mass bias (e.g. redshift and/or mass dependence), the pressure profile, and mass function (Salvati et al. 2018; Ruppin et al. 2019a).

The tension between the cluster number counts from the tSZ and the CMB power spectrum arising from fixing the hydrostatic bias to 0.2 in cosmological analysis has led the scientific community to investigate this discrepancy in more detail. Several studies have been made on the HE mass bias. In this work, a simulated data set of almost 260 clusters from the MUSIC simulations is analysed. We focus on the determination of cluster masses assuming the HE hypothesis, at different redshifts, and compare them with the true cluster mass derived from the simulation data. The HE masses can be estimated using two different equations (Pratt et al. 2019), depending on which ICM thermodynamic quantities are used, i.e. the pressure and electron density (hereafter referred to as SZ mass) or the temperature and electron density (hereafter X-ray mass). There are also different ways of computing the radial gradients of these quantities. In this work, we present a complete study of the HE mass derivations from the different assumptions. Moreover, we also correlate the results for the mass bias with the dynamical state of the MUSIC clusters. A correction to the HE mass, which takes into account the non-thermal pressure contribution arising from gas motions in the ICM, is applied.

We repeated the same analysis for a subsample of 32 objects from MUSIC, in the redshift range 0.5 < *z* < 0.9. This sample, named the NIKA2 twin sample, has been selected to be representative of the clusters observed in the NIKA2 SZ Large Program (LPSZ; Mayet et al. 2020). The LPSZ uniquely exploits the excellent match in sensitivity and spatial resolution of *XMM–Newton* and the NIKA2 camera, which is a millimetre camera installed at the 30-m radio telescope of the Institut de Radioastronomie Millimétrique (IRAM) in Pico Veleta, Spain (Adam et al. 2018; Perotto et al. 2020). A previous analysis of the gas pressure profiles reconstruction on this sample by applying the NIKA2 data reduction pipeline has been already presented in Ruppin et al. (2019b).

The paper is organized as follows. In Section 2, we describe the HE, and the estimate of the cluster mass from different approaches. In Section 3, we review the estimates of the HE mass bias from different numerical hydrodynamical simulations that have been published in literature, showing the relatively large variations in the bias results from the different simulation suites. In Section 4, we briefly introduce the MUSIC simulations and the data used in this analysis, classifying clusters by their dynamical state. The ICM profiles are presented in
SZ effect is a redshift independent probe that consists of the spectral distortion of the cosmic microwave background (CMB) radiation due to the inverse Compton scattering of the CMB photons with the hot electrons of the ICM (for a review, see Carlstrom, Holder & Reese 2002; Kitayama 2014; Mroczkowski et al. 2019). It can be used to directly measure the ICM pressure distribution. The latter can be combined with the electron density from an X-ray observation to infer the cluster mass profile. This method avoids the necessity of obtaining the deep X-ray observations required to measure a spatially resolved temperature profile. To date, the SZ effect induced by thermal electrons (tSZ) has been detected for more than a thousand galaxy clusters, including more than 200 new clusters previously unknown by any other observational means (Kitayama 2014) thanks to new observing facilities such as the South Pole Telescope (SPT; Williamson et al. 2011; Reichardt et al. 2013; Bleem et al. 2020), the Atacama Cosmology Telescope (ACT; Marriage et al. 2010; Hasselfield et al. 2013; Hilton et al. 2018, 2020), and the Planck satellite (Planck Collaboration VIII, XXIX, XXVII 2011a, 2014c, 2016b).

Both X-ray and SZ observations can be used to infer the mass of a cluster. In particular, from the former we can exploit the temperature and the electron density, from the latter the pressure profile. In order to use this information, we need to make three fundamental assumptions: the gas must trace the cluster potential well, it must be spherically symmetric, and in hydrostatic equilibrium (HE; Kravtsov & Borgani 2012 for a review). The mass inferred through this method is called hydrostatic mass (see Section 2). An easy way to track the error made when using HE is the hydrostatic mass bias *b* (Planck Collaboration XX 2014b). It is defined as the difference of the cluster total mass to the one estimated by HE, divided by the total mass *b* = ($M_{tot} - M_{HE}$)/$M_{tot}$.

Galaxy clusters, and, in particular, their number counts are a fundamental cosmological tool. The abundance of clusters and its evolution with redshift are particularly sensitive to the cosmic matter density, $\Omega_m$, and the present amplitude of density fluctuations, characterized by $\sigma_8$, i.e. the rms linear overdensity in spheres of radius $8 \ h^{-1} \ \text{Mpc}$. The CMB primary anisotropies, on the other hand, are related to the density perturbation power spectrum at the time of recombination. A comparison of the amplitude of density perturbations from recombination until today, allows us to look for possible extensions to the concordance $\Lambda$CDM model, such as non-minimal neutrino masses or non-zero curvature contributions (Planck Collaboration XXIV 2016a; Salvati, Douspis & Aghanim 2018).

In Planck Collaboration XX (2014b), a tension between the number of clusters detected by SZ signal, and the number of clusters predicted from the cosmological parameters inferred from the primary CMB power spectrum is reported. This issue was later confirmed in Planck Collaboration XXIV (2016a) and arises when the HE mass bias is fixed to a constant value of 0.2 in cluster cosmological analyses. The hydrostatic mass bias *b* plays an important role in the number counts, because it leads to a modification of the cluster population at a given mass. In particular, it significantly affects the value of $\sigma_8$: the lower the ($1 - b$), the higher the $\sigma_8$, (Salvati et al. 2018; Ruppin et al. 2019a). Salvati et al. (2018) published an update of the constraints on cosmological parameters from the clusters observed by Planck. They find that the bias needed to reconcile CMB constraints with those from the tSZ number counts is ($1 - b$) = 0.62 ± 0.07, which is compatible with the value ($1 - b$) = 0.58 ± 0.04 found in Planck Collaboration XXIV (2016a). This value is confirmed by Koukoufilippas et al. (2020), who cross-correlate Planck maps of the tSZ Compton-y parameter with the galaxy distribution.

These values derived from observations are nevertheless in disagreement with the value of ($1 - b$) estimated from simulations, which is about 0.8. This topic will be deeply discussed later in this work. Moreover, there is a factor of 2.5 more clusters predicted than observed when taking into account the CMB cosmology and a value ($1 - b$) of 0.8 (Andreon 2014; Planck Collaboration XXIV 2016a; Salvati et al. 2018). This value of the bias, derived from simulations, has been recently confirmed by Makiya, Hikage & Komatsu (2020). They perform a joint analysis of power spectra of the tSZ and the cosmic weak lensing shear, in order to obtain a ($1 - b$) constraint which is independent from the primordial CMB spectrum. They find ($1 - b$) = 0.73$^{+0.08}_{-0.13}$, and conclude that the late-time probe (tSZ and cosmic weak lensing shear) cosmologies are consistent with each other, but they could not be totally consistent with the CMB cosmology, which is leading to a different value of the mass bias.

However, the cluster number counts are limited by systematic effects, in particular those affecting the mass estimates. The tSZ power spectrum, in turn, is not measured with sufficient accuracy, especially at small angular scales, to reduce the tension with the CMB. The tSZ cosmological analysis can be improved by considering more realistic and complex hypotheses on the mass bias (e.g. redshift and/or mass dependence), the pressure profile, and mass function (Salvati et al. 2018; Ruppin et al. 2019a).

The tension between the cluster number counts from the tSZ and the CMB power spectrum arising from fixing the hydrostatic bias to 0.2 in cosmological analysis has led the scientific community to investigate this discrepancy in more detail. Several studies have been made on the HE mass bias. In this work, a simulated data set of almost 260 clusters from the MUSIC simulations is analysed. We focus on the determination of cluster masses assuming the HE hypothesis, at different redshifts, and compare them with the true cluster mass derived from the simulation data. The HE masses can be estimated using two different equations (Pratt et al. 2019), depending on which ICM thermodynamic quantities are used, i.e. the pressure and electron density (hereafter referred to as SZ mass) or the temperature and electron density (hereafter X-ray mass). There are also different ways of computing the radial gradients of these quantities. In this work, we present a complete study of the HE mass derivations from the different assumptions. Moreover, we also correlate the results for the mass bias with the dynamical state of the MUSIC clusters. A correction to the HE mass, which takes into account the non-thermal pressure contribution arising from gas motions in the ICM, is applied.

We repeated the same analysis for a subsample of 32 objects from MUSIC, in the redshift range 0.5 < *z* < 0.9. This sample, named the NIKA2 twin sample, has been selected to be representative of the clusters observed in the NIKA2 SZ Large Program (LPSZ; Mayet et al. 2020). The LPSZ uniquely exploits the excellent match in sensitivity and spatial resolution of XMM-Newton and the NIKA2 camera, which is a millimetre camera installed at the 30-m radio telescope of the Institut de Radioastronomie Millimétrique (IRAM) in Pico Veleta, Spain (Adam et al. 2018; Perotto et al. 2020). A previous analysis of the gas pressure profiles reconstruction on this sample by applying the NIKA2 data reduction pipeline has been already presented in Ruppin et al. (2019b).

The paper is organized as follows. In Section 2, we describe the HE, and the estimate of the cluster mass from different approaches. In Section 3, we review the estimates of the HE mass bias from different numerical hydrodynamical simulations that have been published in literature, showing the relatively large variations in the bias results from the different simulation suites. In Section 4, we briefly introduce the MUSIC simulations and the data used in this analysis, classifying clusters by their dynamical state. The ICM profiles are presented in
Section 5. The distribution of the results for HE masses and their biases are described in Section 6, along with the modelling of the non-thermal correction. In Section 7 and Appendix A, we focus on the HE masses and biases for the NIKA2 twin sample. Finally, in Section 8, the main conclusions from these analyses are given.

## 2 HYDROSTATIC EQUILIBRIUM

We use the HE hypothesis to estimate the mass of each clusters (see Kravtsov & Borgani 2012 for a review). It assumes that the gas thermal pressure is balanced by the gravitational force, so that the cluster is in equilibrium. Further assuming that the system is in spherical symmetry and the gas pressure is purely thermal, the total mass inside a sphere of radius r, can be written as

$$M_{\text{HE}}(r) = -\frac{r^2}{G\mu m_p n_e(r)} \frac{dP_{\text{th}}(r)}{dr}, \qquad (1)$$

where G is the gravitational constant, μ is the mean molecular weight of the ICM, here 0.59, m_p is the proton mass, and n_e and P_th are the numerical electron density and the thermal pressure of the gas, respectively. Assuming the equation of state of an ideal gas, it follows that the cluster mass can also be derived from the electron density and temperature T(r) profiles, as

$$M_{\text{HE}}(r) = -\frac{r k_{\text{B}} T(r)}{G \mu m_p} \left[ \frac{d \ln n_e(r)}{d \ln r} + \frac{d \ln T(r)}{d \ln r} \right]. \qquad (2)$$

We refer to equation (1) as M_HE,SZ given that the pressure is estimated by SZ observations, while we refer to equation (2) as M_HE,X because the temperature and the electron density are usually estimated from X-ray observations (Ansarifard et al. 2020).

Deviations from equilibrium could have an impact on observable properties of clusters and may cause systematic errors when equations (1) and (2) are used to estimate cluster masses. A direct comparison with the real mass is usually quantified in the form of a hydrostatic mass bias b (see references in Section 3). It should be stressed that the mass bias can be estimated only if the exact cluster total mass is available, which is the case in simulations. The real mass M_true of a simulated cluster can be easily computed by summing all the dark matter, stars, and gas particle/cell masses inside an aperture radius. The mass bias b_SZ or b_X, at a specific radius, is defined as

$$b = \frac{M_{\text{true}} - M_{\text{HE}}}{M_{\text{true}}}. \qquad (3)$$

The bias defined in this way is usually a positive quantity, since the HE mass often underestimates the true mass. It can happen also in the contrary, leading to a negative bias. Sometimes in literature the opposite difference between the masses is chosen.

## 3 HYDROSTATIC MASS BIAS, STATE OF THE ART

The comparison between tSZ cluster number counts and CMB Planck results has led the community to carefully account for the impact of the HE mass bias on cosmological constraints. Here, we focus on the HE mass bias b, equation (3), computed at R_500. This parameter has been extensively studied with a variety of numerical hydrodynamical simulations.

### 3.1 Previous work in literature

In Fig. 1, we present a compilation of results published in the literature, including the error estimates of the mean values for 1 – b. The mean values are represented as vertical white rectangles, the errors corresponding to each value are the blue coloured regions. The different shades of blue represent the physical processes included in each simulation: non-radiative (light blue), cooling + star formation + supernovae feedback (medium blue), and those including also supermassive black hole (SMBH) feedback (dark blue). In the vertical axis, after the authors reference, we indicate, in parenthesis, the type of bias estimated in each work: SZ for HE masses derived from pressure and density profiles and X-ray for HE masses computed from temperature and density profiles, see Section 2.

Some authors define a positive bias, see the definition in equation (3), others negative, but in Fig. 1, in order to compare them, all biases are taken as positive. Here, we represent the biases from the mass-weighted temperature profiles, in order to focus mainly on the degree of HE bias, since the spectroscopic-like temperature profiles (Mazzotta et al. 2004) are sensitive, in addition, to observational biases in the derived gas density and temperature profiles.

As it can be seen in the Figure, most of the published results for 1 – b are in the range 0.75–0.9, or b ~ 0.25 to 0.10. The majority of these results are in disagreement with respect to the bias needed to reconcile CMB constraints with those from the tSZ number counts of (1 – b) = 0.58 ± 0.04, shown as the vertical red dashed line and shaded region (Planck Collaboration XX, XXIV 2014b, 2016a). Only 3 bias values over 22 are compatible with it within 1σ.

The mean and standard deviation are reported for all the works except for Biffi et al. (2016), where median and MAD (median absolute deviation, computed as median(|b_i – median(b)|)/n-clusters) are presented. We find a X-ray bias of b_X = 0.14+0.11-0.13, using the median and the percentiles (16th and 84th). Also Gupta et al. (2017) report median and percentiles, while Barnes et al. (2020) report the bootstrap errors. There are also some authors, represented without the vertical white rectangle, who give a range of values for b but not a central value with an error (Ameglio et al. 2009; Kay et al. 2012; Martizzi & Agrusa 2016; Henson et al. 2017; Le Brun et al. 2017; Pearce et al. 2020).

### 3.2 Bias dependence on simulation and sample properties

This compilation shows a wide spread in the determination of the bias parameter and their errors. It is obvious that these differences might be attributed to the particularities of the simulations used in each work. To shed some light on this problem, we also account for features like mass resolutions (Fig. 2, central panel); the range of analysed halo masses (Fig. 2, left-hand panel), and the statistics of the total number of clusters (Fig. 2, right-hand panel). Moreover, the physical processes included in each simulation should also be considered. They are represented in Figs 1 and 2 by the different shades of blue of the bars, as already explained above.

In Fig. 2, we indicate the type of code used in each simulation in parenthesis next to the reference of the study. The majority of them are based on different flavours of the smoothed particle hydrodynamics (SPH). Lau, Kravtsov & Nagai (2009), Nelson et al. (2014a), Shi et al. (2016), and Martizzi & Agrusa (2016) use finite volume Eulerian hydrodynamics with adaptive mesh refinement algorithms. The Illustris simulations (Barnes et al. 2020), using the moving mesh code (MM) AREPO (Springel 2010) have the lowest Dark Matter (DM) particle mass, along with the RAMSES code used by Martizzi & Agrusa (2016). Gupta et al. (2017) with Magneticum

<footnote n="1">¹The radius where the cluster density is 500 times the Universe critical density ρ_c at that time, ρ_c = 3H(a)²/(8πG) where H(a) is the Hubble function. M_500 is the mass inside a sphere with radius R_500.</footnote>
Section 5. The distribution of the results for HE masses and their biases are described in Section 6, along with the modelling of the non-thermal correction. In Section 7 and Appendix A, we focus on the HE masses and biases for the NIKA2 twin sample. Finally, in Section 8, the main conclusions from these analyses are given.

## 2 HYDROSTATIC EQUILIBRIUM

We use the HE hypothesis to estimate the mass of each clusters (see Kravtsov & Borgani 2012 for a review). It assumes that the gas thermal pressure is balanced by the gravitational force, so that the cluster is in equilibrium. Further assuming that the system is in spherical symmetry and the gas pressure is purely thermal, the total mass inside a sphere of radius $r$, can be written as

$$M_{\text{HE}}(r) = - \frac{r^2}{G \mu m_p n_e(r)} \frac{dP_{\text{th}}(r)}{dr}, \quad (1)$$

where $G$ is the gravitational constant, $\mu$ is the mean molecular weight of the ICM, here 0.59, $m_p$ is the proton mass, and $n_e$ and $P_{\text{th}}$ are the numerical electron density and the thermal pressure of the gas, respectively. Assuming the equation of state of an ideal gas, it follows that the cluster mass can also be derived from the electron density and temperature $T(r)$ profiles, as

$$M_{\text{HE}}(r) = - \frac{r k_B T(r)}{G \mu m_p} \left[ \frac{d \ln n_e(r)}{d \ln r} + \frac{d \ln T(r)}{d \ln r} \right]. \quad (2)$$

We refer to equation (1) as $M_{\text{HE,SZ}}$ given that the pressure is estimated by SZ observations, while we refer to equation (2) as $M_{\text{HE,X}}$ because the temperature and the electron density are usually estimated from X-ray observations (Ansarifard et al. 2020).

Deviations from equilibrium could have an impact on observable properties of clusters and may cause systematic errors when equations (1) and (2) are used to estimate cluster masses. A direct comparison with the real mass is usually quantified in the form of a hydrostatic mass bias $b$ (see references in Section 3). It should be stressed that the mass bias can be estimated only if the exact cluster total mass is available, which is the case in simulations. The real mass $M_{\text{true}}$ of a simulated cluster can be easily computed by summing all the dark matter, stars, and gas particle/cell masses inside an aperture radius. The mass bias $b_{\text{SZ}}$ or $b_{\text{X}}$, at a specific radius, is defined as

$$b = \frac{M_{\text{true}} - M_{\text{HE}}}{M_{\text{true}}}. \quad (3)$$

The bias defined in this way is usually a positive quantity, since the HE mass often underestimates the true mass. It can happen also in the contrary, leading to a negative bias. Sometimes in literature the opposite difference between the masses is chosen.

## 3 HYDROSTATIC MASS BIAS, STATE OF THE ART

The comparison between tSZ cluster number counts and CMB Planck results has led the community to carefully account for the impact of the HE mass bias on cosmological constraints. Here, we focus on the HE mass bias $b$, equation (3), computed at $R_{500}^1$. This parameter has been extensively studied with a variety of numerical hydrodynamical simulations.

### 3.1 Previous work in literature

In Fig. 1, we present a compilation of results published in the literature, including the error estimates of the mean values for 1 – $b$. The mean values are represented as vertical white rectangles, the errors corresponding to each value are the blue coloured regions. The different shades of blue represent the physical processes included in each simulation: non-radiative (light blue), cooling + star formation + supernovae feedback (medium blue), and those including also supermassive black hole (SMBH) feedback (dark blue). In the vertical axis, after the authors reference, we indicate, in parenthesis, the type of bias estimated in each work: SZ for HE masses derived from pressure and density profiles and X-ray for HE masses computed from temperature and density profiles, see Section 2.

Some authors define a positive bias, see the definition in equation (3), others negative, but in Fig. 1, in order to compare them, all biases are taken as positive. Here, we represent the biases from the mass-weighted temperature profiles, in order to focus mainly on the degree of HE bias, since the spectroscopic-like temperature profiles (Mazzotta et al. 2004) are sensitive, in addition, to observational biases in the derived gas density and temperature profiles.

As it can be seen in the Figure, most of the published results for 1 – $b$ are in the range 0.75–0.9, or $b \sim 0.25$ to 0.10. The majority of these results are in disagreement with respect to the bias needed to reconcile CMB constraints with those from the tSZ number counts of $(1 - b) = 0.58 \pm 0.04$, shown as the vertical red dashed line and shaded region (Planck Collaboration XX, XXIV 2014b, 2016a). Only 3 bias values over 22 are compatible with it within $1\sigma$.

The mean and standard deviation are reported for all the works except for Biffi et al. (2016), where median and MAD (median absolute deviation, computed as median($|b_i - \text{median}(b)|$)/$n_{\text{clusters}}$) are presented. We find a X-ray bias of $b_X = 0.14^{+0.11}_{-0.13}$, using the median and the percentiles (16th and 84th). Also Gupta et al. (2017) report median and percentiles, while Barnes et al. (2020) report the bootstrap errors. There are also some authors, represented without the vertical white rectangle, who give a range of values for $b$ but not a central value with an error (Ameglio et al. 2009; Kay et al. 2012; Martizzi & Agrusa 2016; Henson et al. 2017; Le Brun et al. 2017; Pearce et al. 2020).

### 3.2 Bias dependence on simulation and sample properties

This compilation shows a wide spread in the determination of the bias parameter and their errors. It is obvious that these differences might be attributed to the particularities of the simulations used in each work. To shed some light on this problem, we also account for features like mass resolutions (Fig. 2, central panel); the range of analysed halo masses (Fig. 2, left-hand panel), and the statistics of the total number of clusters (Fig. 2, right-hand panel). Moreover, the physical processes included in each simulation should also be considered. They are represented in Figs 1 and 2 by the different shades of blue of the bars, as already explained above.

In Fig. 2, we indicate the type of code used in each simulation in parenthesis next to the reference of the study. The majority of them are based on different flavours of the smoothed particle hydrodynamics (SPH). Lau, Kravtsov & Nagai (2009), Nelson et al. (2014a), Shi et al. (2016), and Martizzi & Agrusa (2016) use finite volume Eulerian hydrodynamics with adaptive mesh refinement algorithms. The Illustris simulations (Barnes et al. 2020), using the moving mesh code (MM) AREPO (Springel 2010) have the lowest Dark Matter (DM) particle mass, along with the RAMSES code used by Martizzi & Agrusa (2016). Gupta et al. (2017) with Magneticum

<footnote n="1">¹The radius where the cluster density is 500 times the Universe critical density $\rho_c$ at that time, $\rho_c = 3H(a)^2/(8\pi G)$ where $H(a)$ is the Hubble function. $M_{500}$ is the mass inside a sphere with radius $R_{500}$.</footnote>
<figure><img src="image_2.png" /></figure>

**Figure 1.** A compilation from literature of the (1 − *b*) mean values (white rectangle) and their errors (blue bars). They are derived from different hydrodynamical simulation works using the mass-weighted temperature profiles. The vertical red dashed line and its shaded area are the value and the error (1 − *b*) = 0.58 ± 0.04, needed to reconcile the cosmological constraints obtained from *Planck* cluster counts and CMB power spectrum (Planck Collaboration XX, XXIV 2014b, 2016a). The three different shades of blue represent the physics included in each simulations, from the basic set of physical processes, non-radiative in light blue, cooling + star formation + supernovae feedback, in medium blue, to the complete set, including the effects of SMBH feedbacks in dark blue. In the vertical axis, we show the reference for each work and in parenthesis, the type of bias that is derived, from X-Ray, SZ, or both.

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

**Figure 2.** Comparison of the main properties of the simulation works. Left-hand panel: The *M*<sub>500</sub> mass range of clusters used for each simulation work. The colour of the bars have the same meaning as in Fig. 1. Meneghetti et al. (2010) and Martizzi & Agrusa (2016) do not give the exact *M*<sub>500</sub> range, so we give an estimation, from 10<sup>14</sup> to 10<sup>15</sup>. Central panel: The mass resolution of dark matter particles, relative to the MUSIC DM resolution *m*<sub>DM,M</sub> = 1.3 × 10<sup>9</sup> M<sub>☉</sub>, used in this work. Right-hand panel: The number of cluster objects used in each publication, relative to the number of clusters used in this work (*N*<sub>cl,M</sub> = 282). In the vertical axis, next to the reference of each work, we show in parenthesis the computational hydrodynamical method used by each simulation: SPH or AMR. Barnes et al. (2020) employ three different simulations: Illustris (Moving Mesh), BAHAMAS (SPH), and MACSIS (SPH), in this Figure the DM mass from Illustris is represented, since the MACSIS and BAHAMAS codes are also used in Henson et al. (2017) and have the same DM mass.

simulations, and Le Brun et al. (2017) with cosmo-OWLS, employ
the largest number of clusters and have one of the lowest errors on
the bias estimate.

We compared the HE biases in Fig. 1 with the main features of each
simulation shown in Fig. 2. We further plot the biases as a function of

each quantity. For the sake of brevity, we do not show these figures,
but only conclude here that we cannot draw any clear dependence
between the HE bias and the cluster mass range (left-hand panel of
Fig. 2), the particle mass resolution (central panel), or the number of
clusters included in the analysis (right-hand panel).
<figure><img src="image_3.png" /></figure>

**Figure 1.** A compilation from literature of the (1 − *b*) mean values (white rectangle) and their errors (blue bars). They are derived from different hydrodynamical simulation works using the mass-weighted temperature profiles. The vertical red dashed line and its shaded area are the value and the error (1 − *b*) = 0.58 ± 0.04, needed to reconcile the cosmological constraints obtained from *Planck* cluster counts and CMB power spectrum (Planck Collaboration XX, XXIV 2014b, 2016a). The three different shades of blue represent the physics included in each simulations, from the basic set of physical processes, non-radiative in light blue, cooling + star formation + supernovae feedback, in medium blue, to the complete set, including the effects of SMBH feedbacks in dark blue. In the vertical axis, we show the reference for each work and in parenthesis, the type of bias that is derived, from X-Ray, SZ, or both.

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

**Figure 2.** Comparison of the main properties of the simulation works. Left-hand panel: The $M_{500}$ mass range of clusters used for each simulation work. The colour of the bars have the same meaning as in Fig. 1. Meneghetti et al. (2010) and Martizzi & Agrusa (2016) do not give the exact $M_{500}$ range, so we give an estimation, from $10^{14}$ to $10^{15}$. Central panel: The mass resolution of dark matter particles, relative to the MUSIC DM resolution $m_{DM,M} = 1.3 \times 10^9$ M$_\odot$, used in this work. Right-hand panel: The number of cluster objects used in each publication, relative to the number of clusters used in this work ($N_{cl,M} = 282$). In the vertical axis, next to the reference of each work, we show in parenthesis the computational hydrodynamical method used by each simulation: SPH or AMR. Barnes et al. (2020) employ three different simulations: Illustris (Moving Mesh), BAHAMAS (SPH), and MACSIS (SPH), in this Figure the DM mass from Illustris is represented, since the MACSIS and BAHAMAS codes are also used in Henson et al. (2017) and have the same DM mass.

simulations, and Le Brun et al. (2017) with cosmo-OWLS, employ the largest number of clusters and have one of the lowest errors on the bias estimate.

We compared the HE biases in Fig. 1 with the main features of each simulation shown in Fig. 2. We further plot the biases as a function of

each quantity. For the sake of brevity, we do not show these figures, but only conclude here that we cannot draw any clear dependence between the HE bias and the cluster mass range (left-hand panel of Fig. 2), the particle mass resolution (central panel), or the number of clusters included in the analysis (right-hand panel).
The simulation physics comprehends a wide range of processes, but here we gather them in three classes: those which includes the gravitational and non-radiative physics (NR; light blue in Fig. 1); the ‘middle’ set, which includes also radiative processes, like cooling, star formation, and supernovae feedback (CSF, medium blue); and the ‘complete’ set, which, in addition, includes also the feedbacks from supermassive black holes (CSF + SMBH, dark blue). The simulations that have the SMBH feedback seem to have lower errors on the bias, such as Kay et al. (2012), Battaglia et al. (2013), Biffi et al. (2016), McCarthy et al. (2016), Martizzi & Agrusa (2016), Gupta et al. (2017), Pearce et al. (2020), Ansarifard et al. (2020), and Barnes et al. (2020). Even in Meneghetti et al. (2010), the bias error is low mostly for two main reasons: their simulations considered 1/3 of Spitzer thermal conductivity that homogenize the medium and the statistics are limited due to the small sample size (see Fig. 2). It is not straightforward to derive the impact of SMBH on the HE mass bias from the different simulation results compiled in Fig. 1, because each simulation has different features, see Fig. 2. In this paper, we are comparing an homogeneous and statistical significant set of clusters simulated with different physics flavours. Thus, we can provide a definitive answer on this issue (see Section 6.3.3).

### 3.3 Bias dependence on other quantities or measurements

The more compelling questions about the HE mass bias are its dependence on redshift, on cluster dynamical state, or on whether the bias is calculated from spectroscopic-like temperature. The redshift dependence has been studied by several authors (Piffaretti & Valdarnini 2008; Lau et al. 2009; Henson et al. 2017; Le Brun et al. 2017). The expectation is that clusters, going at higher redshifts and being less relaxed, tend to have more mass in substructures (Neto et al. 2007; Vazza et al. 2018) risking to violate the HE hypothesis. However, this behaviour has not yet been confirmed in any work.

According to their dynamical state, clusters are usually classified in two main classes: the relaxed ones, well described by spherical symmetry and HE, and disturbed clusters, the opposite. Non-thermal gas motions in the ICM and the non-spherical symmetry of a cluster (usually disturbed) will most likely lead to a larger HE bias and a larger scatter, with respect to the more regular (relaxed) clusters. Several authors analysed the mass bias dependence with the dynamical state, like Piffaretti & Valdarnini (2008), Rasia et al. (2012), Nelson et al. (2014a), Henson et al. (2017), and Ansarifard et al. (2020). They all find no significant distinction between the mass bias of regular and disturbed clusters, given the large dispersion (Cialone et al. 2018). However, Biffi et al. (2016) observe that cool-core (CC) and non-cool-core (NCC) clusters behave differently, with a larger bias for NCC, especially in the innermost cluster regions.

The spectroscopic-like temperature, see Section 5.3, is estimated by weighting the temperature by the X-ray emission of each gas particle. Usually the HE mass bias estimated from this temperature by combining it with the electron density is larger than the one from the mass-weighted temperature and shows a dependence with the true mass of the haloes. This is mainly due to temperature inhomogeneities (Rasia et al. 2006, 2012, 2014; Henson et al. 2017; Le Brun et al. 2017; Barnes et al. 2020; Pearce et al. 2020). Moreover, Piffaretti & Valdarnini (2008) find that the spectroscopic bias also depends on the dynamical state of clusters, with large biases found in the most disturbed clusters (see also Biffi et al. 2014).

In the last years, the HE mass biases (SZ and X-ray) were often studied together (Ansarifard et al. 2020; Pearce et al. 2020). Ansarifard et al. (2020), for example, analyse more than 300 simulated massive clusters, from ‘The Three Hundred Project’ (Cui et al. 2018).

They find that a robust correction to the HE mass bias can be inferred when the gas inhomogeneity from X-ray maps are combined with the asymptotic external slope of the gas density or pressure profiles, which can be derived from X-ray and SZ effect observations. Both SZ and X-ray biases are estimated, with values of 10 per cent, by using models to fit ICM radial profiles.

## 4 THE SIMULATED DATA SET

### 4.1 MUSIC simulations

The clusters analysed in this work are taken from the MUSIC project (Sembolini et al. 2013) that consists of two sets of resimulated clusters extracted from two large volume simulations: the MUSIC-1 sample, extracted from the 500 $h^{-1}$ Mpc MareNostrum Universe simulation box (Gottlöber & Yepes 2007), and the MUSIC-2 sample, from the $1 h^{-1}$ Gpc MultiDark (MD) simulation box (Prada et al. 2012).

In this work, the 258 zoomed regions around the most massive clusters in the MUSIC-2 data base were analysed. From a low-resolution version (256³ particles) of the two simulations, the particles inside a sphere of 6 $h^{-1}$ Mpc radius at $z=0$ are mapped back to the initial redshift, using the Klypin et al. (2001) zooming technique, to identify their corresponding Lagrangian regions. These regions are then resimulated with high resolution and populated with SPH gas particles. The original MD dark-matter-only simulation was performed with L-GADGET2 code (Klypin et al. 2016) and adopting WMAP7 + BAO + SNI cosmology: $\Omega_m = 0.27$, $\Omega_b = 0.0469$, $\Omega_\Lambda = 0.73$, $\sigma_8 = 0.82$, $n = 0.95$, and $h = 0.7$ (Komatsu et al. 2011).

All the resimulations are done with the TreePM + SPH GADGET code (Springel 2005), and include three different classes of physical processes, labelled as *flavours*: non-radiative (NR), cooling and star formation (CSF), and active galactic nuclei (AGNs). The NR flavour includes only the gravitational and gasdynamical effects, while the CSF flavour includes radiative processes, like star formation, feedback from supernovae, UV photoionization, and radiative cooling (Sembolini et al. 2013), and finally AGNs, where the AGNs and their feedback are added, using the models for SMBH feedbacks (Planelles et al. 2013).

The clusters were identified using a Bound Density Maxima halo finder (see also AHF halo finder; Knollmann & Knebe 2009). Since we take all the objects above a given mass, the MUSIC-2 catalogue constitutes a complete volume limited sample. Our cluster masses $M_{500}$ range from $1.9 \times 10^{14} M_{\odot}$ to $1.7 \times 10^{15} M_{\odot}$ at $z=0$. All of them were resimulated with NR, CSF, and AGN flavours with a DM and gas mass resolution of $m_{\text{DM}} = 1.29 \times 10^9 M_{\odot}$ and $m_{\text{gas}} = 2.7 \times 10^8 M_{\odot}$. MUSIC-2 cluster regions have been saved at specific redshifts, so in this analysis we also studied the cosmic evolution of these objects, using only the main progenitors of the $z=0$ clusters at redshifts 0.11, 0.33, 0.43, 0.54, 0.67, and 0.82.

### 4.2 The NIKA2 LPSZ twin sample

From the MUSIC-2 data set, a subsample of objects has been selected to reproduce the clusters present in the NIKA2 Large Program SZ (LPSZ) catalogue (Mayet et al. 2020). NIKA2 (Calvo et al. 2016; Adam et al. 2018; Perotto et al. 2020) is the new multipixels camera at 150 and 260 GHz installed at the 30-m telescope of the Institut de Radioastronomie Millimétrique (IRAM). The NIKA2 SZ large program consists of mapping the tSZ signal of a representative sample of 50 galaxy clusters at high angular resolution (18 and 11 arcsec in two bands) and in the $0.5 < z < 0.9$ redshift range. The cluster sample
The simulation physics comprehends a wide range of processes, but here we gather them in three classes: those which includes the gravitational and non-radiative physics (NR; light blue in Fig. 1); the ‘middle’ set, which includes also radiative processes, like cooling, star formation, and supernovae feedback (CSF, medium blue); and the ‘complete’ set, which, in addition, includes also the feedbacks from supermassive black holes (CSF + SMBH, dark blue). The simulations that have the SMBH feedback seem to have lower errors on the bias, such as Kay et al. (2012), Battaglia et al. (2013), Biffi et al. (2016), McCarthy et al. (2016), Martizzi & Agrusa (2016), Gupta et al. (2017), Pearce et al. (2020), Ansarifard et al. (2020), and Barnes et al. (2020). Even in Meneghetti et al. (2010), the bias error is low mostly for two main reasons: their simulations considered 1/3 of Spitzer thermal conductivity that homogenize the medium and the statistics are limited due to the small sample size (see Fig. 2). It is not straightforward to derive the impact of SMBH on the HE mass bias from the different simulation results compiled in Fig. 1, because each simulation has different features, see Fig. 2. In this paper, we are comparing an homogeneous and statistical significant set of clusters simulated with different physics flavours. Thus, we can provide a definitive answer on this issue (see Section 6.3.3).

## 3.3 Bias dependence on other quantities or measurements

The more compelling questions about the HE mass bias are its dependence on redshift, on cluster dynamical state, or on whether the bias is calculated from spectroscopic-like temperature. The redshift dependence has been studied by several authors (Piffaretti & Valdarnini 2008; Lau et al. 2009; Henson et al. 2017; Le Brun et al. 2017). The expectation is that clusters, going at higher redshifts and being less relaxed, tend to have more mass in substructures (Neto et al. 2007; Vazza et al. 2018) risking to violate the HE hypothesis. However, this behaviour has not yet been confirmed in any work.

According to their dynamical state, clusters are usually classified in two main classes: the relaxed ones, well described by spherical symmetry and HE, and disturbed clusters, the opposite. Non-thermal gas motions in the ICM and the non-spherical symmetry of a cluster (usually disturbed) will most likely lead to a larger HE bias and a larger scatter, with respect to the more regular (relaxed) clusters. Several authors analysed the mass bias dependence with the dynamical state, like Piffaretti & Valdarnini (2008), Rasia et al. (2012), Nelson et al. (2014a), Henson et al. (2017), and Ansarifard et al. (2020). They all find no significant distinction between the mass bias of regular and disturbed clusters, given the large dispersion (Cialone et al. 2018). However, Biffi et al. (2016) observe that cool-core (CC) and non-cool-core (NCC) clusters behave differently, with a larger bias for NCC, especially in the innermost cluster regions.

The spectroscopic-like temperature, see Section 5.3, is estimated by weighting the temperature by the X-ray emission of each gas particle. Usually the HE mass bias estimated from this temperature by combining it with the electron density is larger than the one from the mass-weighted temperature and shows a dependence with the true mass of the haloes. This is mainly due to temperature inhomogeneities (Rasia et al. 2006, 2012, 2014; Henson et al. 2017; Le Brun et al. 2017; Barnes et al. 2020; Pearce et al. 2020). Moreover, Piffaretti & Valdarnini (2008) find that the spectroscopic bias also depends on the dynamical state of clusters, with large biases found in the most disturbed clusters (see also Biffi et al. 2014).

In the last years, the HE mass biases (SZ and X-ray) were often studied together (Ansarifard et al. 2020; Pearce et al. 2020). Ansarifard et al. (2020), for example, analyse more than 300 simulated massive clusters, from ‘The Three Hundred Project’ (Cui et al. 2018).

They find that a robust correction to the HE mass bias can be inferred when the gas inhomogeneity from X-ray maps are combined with the asymptotic external slope of the gas density or pressure profiles, which can be derived from X-ray and SZ effect observations. Both SZ and X-ray biases are estimated, with values of 10 per cent, by using models to fit ICM radial profiles.

# 4 THE SIMULATED DATA SET

## 4.1 MUSIC simulations

The clusters analysed in this work are taken from the MUSIC project (Sembolini et al. 2013) that consists of two sets of resimulated clusters extracted from two large volume simulations: the MUSIC-1 sample, extracted from the 500 $h^{-1}$ Mpc MareNostrum Universe simulation box (Gottlöber & Yepes 2007), and the MUSIC-2 sample, from the $1 h^{-1}$ Gpc MultiDark (MD) simulation box (Prada et al. 2012).

In this work, the 258 zoomed regions around the most massive clusters in the MUSIC-2 data base were analysed. From a low-resolution version (256³ particles) of the two simulations, the particles inside a sphere of 6 $h^{-1}$ Mpc radius at $z=0$ are mapped back to the initial redshift, using the Klypin et al. (2001) zooming technique, to identify their corresponding Lagrangian regions. These regions are then resimulated with high resolution and populated with SPH gas particles. The original MD dark-matter-only simulation was performed with L-GADGET2 code (Klypin et al. 2016) and adopting WMAP7 + BAO + SNI cosmology: $\Omega_m = 0.27$, $\Omega_b = 0.0469$, $\Omega_\Lambda = 0.73$, $\sigma_8 = 0.82$, $n = 0.95$, and $h = 0.7$ (Komatsu et al. 2011).

All the resimulations are done with the TreePM + SPH GADGET code (Springel 2005), and include three different classes of physical processes, labelled as *flavours*: non-radiative (NR), cooling and star formation (CSF), and active galactic nuclei (AGNs). The NR flavour includes only the gravitational and gasdynamical effects, while the CSF flavour includes radiative processes, like star formation, feedback from supernovae, UV photoionization, and radiative cooling (Sembolini et al. 2013), and finally AGNs, where the AGNs and their feedback are added, using the models for SMBH feedbacks (Planelles et al. 2013).

The clusters were identified using a Bound Density Maxima halo finder (see also AHF halo finder; Knollmann & Knebe 2009). Since we take all the objects above a given mass, the MUSIC-2 catalogue constitutes a complete volume limited sample. Our cluster masses $M_{500}$ range from $1.9 \times 10^{14} M_{\odot}$ to $1.7 \times 10^{15} M_{\odot}$ at $z=0$. All of them were resimulated with NR, CSF, and AGN flavours with a DM and gas mass resolution of $m_{\text{DM}} = 1.29 \times 10^9 M_{\odot}$ and $m_{\text{gas}} = 2.7 \times 10^8 M_{\odot}$. MUSIC-2 cluster regions have been saved at specific redshifts, so in this analysis we also studied the cosmic evolution of these objects, using only the main progenitors of the $z=0$ clusters at redshifts 0.11, 0.33, 0.43, 0.54, 0.67, and 0.82.

## 4.2 The NIKA2 LPSZ twin sample

From the MUSIC-2 data set, a subsample of objects has been selected to reproduce the clusters present in the NIKA2 Large Program SZ (LPSZ) catalogue (Mayet et al. 2020). NIKA2 (Calvo et al. 2016; Adam et al. 2018; Perotto et al. 2020) is the new multipixels camera at 150 and 260 GHz installed at the 30-m telescope of the Institut de Radioastronomie Millimétrique (IRAM). The NIKA2 SZ large program consists of mapping the tSZ signal of a representative sample of 50 galaxy clusters at high angular resolution (18 and 11 arcsec in two bands) and in the $0.5 < z < 0.9$ redshift range. The cluster sample
<table><caption>**Table 1.** Fractions of relaxed, disturbed, and intermediate clusters at each redshift for the MUSIC sample and for the NIKA2 subsample.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">MUSIC</th><th colspan="3">NIKA2</th></tr><tr><th>Rel. (per cent)</th><th>Interm. (per cent)</th><th>Dist. (per cent)</th><th>Rel. (per cent)</th><th>Interm. (per cent)</th><th>Dist. (per cent)</th></tr></thead><tbody><tr><td>0.00</td><td>52</td><td>29</td><td>19</td><td></td><td></td><td></td></tr><tr><td>0.11</td><td>50</td><td>34</td><td>16</td><td></td><td></td><td></td></tr><tr><td>0.33</td><td>53</td><td>28</td><td>19</td><td></td><td></td><td></td></tr><tr><td>0.43</td><td>49</td><td>34</td><td>17</td><td></td><td></td><td></td></tr><tr><td>0.54</td><td>47</td><td>30</td><td>23</td><td>44</td><td>39</td><td>17</td></tr><tr><td>0.67</td><td>48</td><td>30</td><td>22</td><td></td><td></td><td></td></tr><tr><td>0.82</td><td>52</td><td>30</td><td>18</td><td>57</td><td>22</td><td>21</td></tr></tbody></table>

was extracted from the tSZ catalogues established by the Planck and ACT collaborations (Hasselfield et al. 2013; Planck Collaboration XXVII 2016b), and the selected clusters homogeneously populate the mass range with $M_{500} > 3 \times 10^{14} M_{\odot}$ (Mayet et al. 2020) and redshift range. The MUSIC NIKA2 twin sample closely matches the same mass-redshift space as the NIKA2 tSZ large program. For the redshift bin $0.5 < z < 0.7$, 18 clusters were chosen from the MUSIC-2 catalogue at redshift 0.54. For the $0.7 < z < 0.9$ bin, 14 clusters from MUSIC redshift 0.82 were also selected. The same mass cut applied to the Planck and ACT catalogue of $M_{500} > 3 \times 10^{14} M_{\odot}$ was applied to the MUSIC sample, so only the clusters with an HE mass estimate above this threshold are used, in order to make the two sample comparable. A previous analysis on the gas pressure profiles recovered from the NIKA2 twin sample has been already performed (Ruppin et al. 2019b; De Petris et al. 2020).

## 4.3 Characterization of the cluster dynamical state and morphology

Throughout this analysis, the dynamical state of the clusters has been inferred by two 3D indicators (Neto et al. 2007; Sembolini et al. 2013; Cialone et al. 2018; for a novel approach to infer the cluster morphology and dynamical state with the Zernike polynomials, see Capalbo et al. 2020). Here, we focus on the dynamical state of clusters inside $R_{500}$, which is more in agreement with what is measured in observations. The considered 3D dynamical state estimators are

(i) $M_{sub}/M_{500}$, the ratio between the mass of the most massive cluster substructure and the total cluster mass inside $R_{500}$. This indicator is mainly sensible to strong mergers. Therefore, it is useful to find the really disturbed clusters. An alternative definition of $M_{sub}$ would be to account for the mass of all the substructures within the aperture radius. This approach is more sensible to a cluster relaxation state (Cui et al. 2017; De Luca et al. 2020).

(ii) $\Delta_r$, the offset between the central density peak, $r_\delta$, and the centre of mass of the cluster, $r_{cm}$, normalized to the aperture radius $R_{500}$:

$$ \Delta_r = \frac{|\mathbf{r}_\delta - \mathbf{r}_{cm}|}{R_{500}}. \qquad (4) $$

In order to have a relaxed cluster, both indicators should be smaller than a given threshold, which varies depending on the authors (D'Onghia & Navarro 2007; Macciò et al. 2007). Here, following Cialone et al. (2018), both indicators should be lower than 0.1 to define a cluster as relaxed, and greater than 0.1 to have a disturbed one. In the cases in which the two indicators provide contradictory answers, the cluster is classified as intermediate or hybrid. According to this selection for the dynamical state estimators, as seen in Table 1, the MUSIC data set contains roughly 50 per cent of relaxed clusters

for all the redshift intervals considered. The same happens for the NIKA2 sub sample.

We can combine the information provided by the two dynamical indicators in a single 'relaxation' parameter $\chi_{DS}$, defined as in Haggar et al. (2020), but dropping the virial ratio $\eta$, in order to keep the same number of relaxed clusters from the definition of the two separate dynamical indicators, $M_{sub}/M_{500}$ and $\Delta_r$ (see De Luca et al. 2020). $\chi_{DS}$ is a continuous, non-binary, estimate of the dynamical state

$$ \chi_{DS} = \sqrt{\frac{2}{\left(\frac{\Delta_r}{0.1}\right)^2 + \left(\frac{M_{sub}/M_{500}}{0.1}\right)^2}}. \qquad (5) $$

All clusters that have this parameter higher than 1 are dynamically relaxed. We will study how this parameter is correlated with $b$ (in Section 6.3.1).

# 5 ICM PROFILES OF THE MUSIC CLUSTERS

For each MUSIC cluster, we compute the 3D radial profiles of the ICM thermodynamic properties. The cluster is divided into spherical shells, from 0 (the core) to $3R_{vir}$ (the outskirts), each with a thickness of 10 kpc. The gas pressure and the electron density are evaluated as the median of all the SPH gas particles inside each spherical shell, while the temperature used is the mass-weighted one. The associated uncertainty of the median value is estimated by MAD

$$ \text{MAD} = \text{median}(|X_i - \text{median}(X)|). \qquad (6) $$

The uncertainties associated with the median profiles usually increase in the cluster outskirts due to the deviations from the spherical symmetry and from a homogenous distribution of the ICM density and temperature (e.g. presence of clumps or disturbances generated by accreting material).

## 5.1 Pressure profile

The cluster pressure profile is well modelled by the general-ized Navarro-Frenk-White (gNFW) model, introduced by Nagai, Kravtsov & Vikhlinin (2007)

$$ \frac{P(r)}{P_{500}} = \frac{P_0}{x^c(1 + x^a)^{\frac{b-c}{a}}} \qquad (7) $$

with $x = r/r_s$ a dimensionless radial distance normalized to the scale radius $r_s = R_{500}/c_{500}$, where $c_{500}$ is the concentration parameter. The parameters $b$ and $c$ are the slopes for outer and inner region radii, respectively, and $a$ is the steepness of the transition between the two regimes. The normalization $P_{500}$ is inferred by the scaling relation between the pressure content and the cluster total mass in the self-similar model (Arnaud et al. 2010, hereafter A10) purely based on gravitation:

$$ P_{500} = 1.65 \times 10^{-3} E_z^{8/3} \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{2/3} h_{70}^2 \text{keV/cm}^3, \qquad (8) $$

$E_z$ is the ratio of the Hubble constant at redshift $z$ to its present value $H_0$, and $h_{70} = H_0/[70 \text{ km s}^{-1} \text{ Mpc}^{-1}]$.

A worthwhile step is to compare our simulated cluster sample with observations. A10 computed the mean pressure profile of galaxy clusters using observed clusters from REXCESS, a representative sample of 33 local clusters ($z < 0.2$) drawn from the REFLEX catalogue and observed with XMM-Newton satellite, and three large samples of simulated clusters at redshift zero extracted from hydrodynamical simulations (Borgani et al. 2004; Nagai et al. 2007;
<table><caption>**Table 1.** Fractions of relaxed, disturbed, and intermediate clusters at each redshift for the MUSIC sample and for the NIKA2 subsample.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">MUSIC</th><th colspan="3">NIKA2</th></tr><tr><th>Rel. (per cent)</th><th>Interm. (per cent)</th><th>Dist. (per cent)</th><th>Rel. (per cent)</th><th>Interm. (per cent)</th><th>Dist. (per cent)</th></tr></thead><tbody><tr><td>0.00</td><td>52</td><td>29</td><td>19</td><td></td><td></td><td></td></tr><tr><td>0.11</td><td>50</td><td>34</td><td>16</td><td></td><td></td><td></td></tr><tr><td>0.33</td><td>53</td><td>28</td><td>19</td><td></td><td></td><td></td></tr><tr><td>0.43</td><td>49</td><td>34</td><td>17</td><td></td><td></td><td></td></tr><tr><td>0.54</td><td>47</td><td>30</td><td>23</td><td>44</td><td>39</td><td>17</td></tr><tr><td>0.67</td><td>48</td><td>30</td><td>22</td><td></td><td></td><td></td></tr><tr><td>0.82</td><td>52</td><td>30</td><td>18</td><td>57</td><td>22</td><td>21</td></tr></tbody></table>

was extracted from the tSZ catalogues established by the Planck and ACT collaborations (Hasselfield et al. 2013; Planck Collaboration XXVII 2016b), and the selected clusters homogeneously populate the mass range with $M_{500} > 3 \times 10^{14} M_{\odot}$ (Mayet et al. 2020) and redshift range. The MUSIC NIKA2 twin sample closely matches the same mass-redshift space as the NIKA2 tSZ large program. For the redshift bin $0.5 < z < 0.7$, 18 clusters were chosen from the MUSIC-2 catalogue at redshift 0.54. For the $0.7 < z < 0.9$ bin, 14 clusters from MUSIC redshift 0.82 were also selected. The same mass cut applied to the Planck and ACT catalogue of $M_{500} > 3 \times 10^{14} M_{\odot}$ was applied to the MUSIC sample, so only the clusters with an HE mass estimate above this threshold are used, in order to make the two sample comparable. A previous analysis on the gas pressure profiles recovered from the NIKA2 twin sample has been already performed (Ruppin et al. 2019b; De Petris et al. 2020).

## 4.3 Characterization of the cluster dynamical state and morphology

Throughout this analysis, the dynamical state of the clusters has been inferred by two 3D indicators (Neto et al. 2007; Sembolini et al. 2013; Cialone et al. 2018; for a novel approach to infer the cluster morphology and dynamical state with the Zernike polynomials, see Capalbo et al. 2020). Here, we focus on the dynamical state of clusters inside $R_{500}$, which is more in agreement with what is measured in observations. The considered 3D dynamical state estimators are

(i) $M_{sub}/M_{500}$, the ratio between the mass of the most massive cluster substructure and the total cluster mass inside $R_{500}$. This indicator is mainly sensible to strong mergers. Therefore, it is useful to find the really disturbed clusters. An alternative definition of $M_{sub}$ would be to account for the mass of all the substructures within the aperture radius. This approach is more sensible to a cluster relaxation state (Cui et al. 2017; De Luca et al. 2020).

(ii) $\Delta_r$, the offset between the central density peak, $r_\delta$, and the centre of mass of the cluster, $r_{cm}$, normalized to the aperture radius $R_{500}$:

$$ \Delta_r = \frac{|\mathbf{r}_\delta - \mathbf{r}_{cm}|}{R_{500}}. \qquad (4) $$

In order to have a relaxed cluster, both indicators should be smaller than a given threshold, which varies depending on the authors (D'Onghia & Navarro 2007; Macciò et al. 2007). Here, following Cialone et al. (2018), both indicators should be lower than 0.1 to define a cluster as relaxed, and greater than 0.1 to have a disturbed one. In the cases in which the two indicators provide contradictory answers, the cluster is classified as intermediate or hybrid. According to this selection for the dynamical state estimators, as seen in Table 1, the MUSIC data set contains roughly 50 per cent of relaxed clusters

for all the redshift intervals considered. The same happens for the NIKA2 sub sample.

We can combine the information provided by the two dynamical indicators in a single 'relaxation' parameter $\chi_{DS}$, defined as in Haggar et al. (2020), but dropping the virial ratio $\eta$, in order to keep the same number of relaxed clusters from the definition of the two separate dynamical indicators, $M_{sub}/M_{500}$ and $\Delta_r$ (see De Luca et al. 2020). $\chi_{DS}$ is a continuous, non-binary, estimate of the dynamical state

$$ \chi_{DS} = \sqrt{\frac{2}{\left(\frac{\Delta_r}{0.1}\right)^2 + \left(\frac{M_{sub}/M_{500}}{0.1}\right)^2}}. \qquad (5) $$

All clusters that have this parameter higher than 1 are dynamically relaxed. We will study how this parameter is correlated with $b$ (in Section 6.3.1).

# 5 ICM PROFILES OF THE MUSIC CLUSTERS

For each MUSIC cluster, we compute the 3D radial profiles of the ICM thermodynamic properties. The cluster is divided into spherical shells, from 0 (the core) to $3R_{vir}$ (the outskirts), each with a thickness of 10 kpc. The gas pressure and the electron density are evaluated as the median of all the SPH gas particles inside each spherical shell, while the temperature used is the mass-weighted one. The associated uncertainty of the median value is estimated by MAD

$$ \text{MAD} = \text{median}(|X_i - \text{median}(X)|). \qquad (6) $$

The uncertainties associated with the median profiles usually increase in the cluster outskirts due to the deviations from the spherical symmetry and from a homogenous distribution of the ICM density and temperature (e.g. presence of clumps or disturbances generated by accreting material).

## 5.1 Pressure profile

The cluster pressure profile is well modelled by the general-ized Navarro–Frenk–White (gNFW) model, introduced by Nagai, Kravtsov & Vikhlinin (2007)

$$ \frac{P(r)}{P_{500}} = \frac{P_0}{x^c(1 + x^a)^{\frac{b-c}{a}}} \qquad (7) $$

with $x = r/r_s$ a dimensionless radial distance normalized to the scale radius $r_s = R_{500}/c_{500}$, where $c_{500}$ is the concentration parameter. The parameters $b$ and $c$ are the slopes for outer and inner region radii, respectively, and $a$ is the steepness of the transition between the two regimes. The normalization $P_{500}$ is inferred by the scaling relation between the pressure content and the cluster total mass in the self-similar model (Arnaud et al. 2010, hereafter A10) purely based on gravitation:

$$ P_{500} = 1.65 \times 10^{-3} E_z^{8/3} \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{2/3} h_{70}^2 \text{keV/cm}^3, \qquad (8) $$

$E_z$ is the ratio of the Hubble constant at redshift $z$ to its present value $H_0$, and $h_{70} = H_0/[70 \text{ km s}^{-1} \text{Mpc}^{-1}]$.

A worthwhile step is to compare our simulated cluster sample with observations. A10 computed the mean pressure profile of galaxy clusters using observed clusters from REXCESS, a representative sample of 33 local clusters ($z < 0.2$) drawn from the REFLEX catalogue and observed with XMM–Newton satellite, and three large samples of simulated clusters at redshift zero extracted from hydrodynamical simulations (Borgani et al. 2004; Nagai et al. 2007;
<table><caption>**Table 2.** The best-fitting parameters and their errors for the gNFW overall median pressure profile at low redshifts ($z < 0.5$). In the first part of the table, the parameters are listed for the three simulation flavours (AGN, CSF, and NR), followed by the parameters from the A10 universal profile and the P13. In the second part of the table, we show the best-fitting parameters for relaxed MUSIC clusters (AGN flavour) and for the CC clusters from the A10 and P13. The description is exactly the same for the third part of the table, but in this case for the disturbed clusters and the NCC clusters from P13.</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>AGN</td><td>8.87 ± 0.52</td><td>1.77 ± 0.04</td><td>1.15 ± 0.03</td><td>4.29 ± 0.06</td><td>0.29 ± 0.02</td></tr><tr><td>CSF</td><td>2.55 ± 0.10</td><td>1.20 ± 0.05</td><td>1.14 ± 0.03</td><td>4.50 ± 0.09</td><td>0.74 ± 0.02</td></tr><tr><td>NR</td><td>20.00 ± 0.20</td><td>2.67 ± 0.13</td><td>1.07 ± 0.05</td><td>3.40 ± 0.09</td><td>0.23 ± 0.03</td></tr><tr><td>P13</td><td>6.41</td><td>1.81</td><td>1.33</td><td>4.13</td><td>0.31</td></tr><tr><td>A10</td><td>8.403</td><td>1.177</td><td>1.051</td><td>5.491</td><td>0.308</td></tr><tr><td>AGN relaxed</td><td>13.40 ± 0.72</td><td>2.39 ± 0.03</td><td>1.24 ± 0.03</td><td>3.95 ± 0.04</td><td>0.27 ± 0.02</td></tr><tr><td>P13 CC</td><td>11.82</td><td>0.60</td><td>0.76</td><td>6.58</td><td>0.31</td></tr><tr><td>A10 CC</td><td>3.25</td><td>1.13</td><td>1.22</td><td>5.49</td><td>0.77</td></tr><tr><td>AGN disturbed</td><td>3.29 ± 0.24</td><td>1.35 ± 0.09</td><td>1.38 ± 0.07</td><td>4.35 ± 0.16</td><td>0.39 ± 0.03]</td></tr><tr><td>P13 NCC</td><td>4.72</td><td>2.19</td><td>1.82</td><td>3.62</td><td>0.31</td></tr><tr><td>A10 dist.</td><td>3.20</td><td>1.08</td><td>1.41</td><td>5.49</td><td>0.38</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The median pressure profile of the MUSIC clusters (almost 1050) at redshift $z < 0.5$ for the AGN (left-hand panel), CSF (central panel), and NR (right-hand panel) simulation flavours. The grey lines represent the individual profiles for each cluster. The overall median profile is represented by black dots, with error bars computed by the MAD (equation 6). The yellow area shows the variation of the best fit for the gNFW between $1\sigma$ and $3\sigma$ intervals. The bottom panel of each plot shows the relative difference between the fit of the MUSIC profile and the corresponding profiles from A10 (green) and P13 (magenta), respectively.</figcaption></figure>

Piffaretti & Valdarnini 2008). The fit on the profile was performed in the radial range [0.03–4] × $R_{500}$. In this radial range, the observed profile is limited to $R_{500}$ and the region outside this radius was extrapolated according to the predictions from numerical simulations. They describe the resulting pressure profiles as universal, since it fits well both data from simulations and observations, with parameters listed in Table 2, fifth row.

The Planck Collaboration (Planck Collaboration V 2014a, hereafter P13) compared the median gNFW pressure profile of 62 nearby ($z < 0.5$) massive observed clusters with the A10 profile, finding that there is a very good agreement in the cluster intermediate radii between the two. However, within the core, i.e. $R < 0.15R_{500}$, the observed profile lies significantly below the A10 profile. The fit was done in the radial range [0.02–3] × $R_{500}$ and the parameters are reported in the fourth row of Table 2.

To compare our simulated sample with observations, we compare the median pressure profiles (dots in Fig. 3) for all the MUSIC cluster sets at low redshifts ($z < 0.5$, a sample of almost 1050 clusters), with the A10 and P13 gNFW profiles, represented in green and violet, respectively. The three panels show the MUSIC

median profiles for the simulation flavours, AGN, CSF, and NR from left to right, with the MAD, see equation (6), as associated uncertainty. The yellow line and the shaded regions represent the MUSIC pressure profile fit, with the gNFW model parameters for AGN, CSF, and NR listed in the first three rows of Table 2. In the bottom plot of each panel, we present the relative difference between the MUSIC pressure profile fit, $f_M$, and the profile from A10 or P13

$$ \Delta = \frac{f_M - f_{A10/P13}}{f_M} \qquad (9) $$

The AGN flavour is the data set showing better agreement with both A10 and P13 profiles, especially for the A10’s universal profile, a very good approximation until $R_{500}$. In the case of the CSF flavour, the MUSIC profiles are steeper starting from $0.1R_{500}$, while in the NR case the MUSIC profiles are higher than the observed profiles within even a larger region (approximately $0.3R_{500}$). Since only the AGN set provides a reliable description of the observed profiles, we used this set to check the redshift dependence of the universal profile, extended to the high-redshift regime ($0.54 < z < 0.82$). As
<table><caption>**Table 2.** The best-fitting parameters and their errors for the gNFW overall median pressure profile at low redshifts ($z < 0.5$). In the first part of the table, the parameters are listed for the three simulation flavours (AGN, CSF, and NR), followed by the parameters from the A10 universal profile and the P13. In the second part of the table, we show the best-fitting parameters for relaxed MUSIC clusters (AGN flavour) and for the CC clusters from the A10 and P13. The description is exactly the same for the third part of the table, but in this case for the disturbed clusters and the NCC clusters from P13.</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>AGN</td><td>8.87 ± 0.52</td><td>1.77 ± 0.04</td><td>1.15 ± 0.03</td><td>4.29 ± 0.06</td><td>0.29 ± 0.02</td></tr><tr><td>CSF</td><td>2.55 ± 0.10</td><td>1.20 ± 0.05</td><td>1.14 ± 0.03</td><td>4.50 ± 0.09</td><td>0.74 ± 0.02</td></tr><tr><td>NR</td><td>20.00 ± 0.20</td><td>2.67 ± 0.13</td><td>1.07 ± 0.05</td><td>3.40 ± 0.09</td><td>0.23 ± 0.03</td></tr><tr><td>P13</td><td>6.41</td><td>1.81</td><td>1.33</td><td>4.13</td><td>0.31</td></tr><tr><td>A10</td><td>8.403</td><td>1.177</td><td>1.051</td><td>5.491</td><td>0.308</td></tr><tr><td>AGN relaxed</td><td>13.40 ± 0.72</td><td>2.39 ± 0.03</td><td>1.24 ± 0.03</td><td>3.95 ± 0.04</td><td>0.27 ± 0.02</td></tr><tr><td>P13 CC</td><td>11.82</td><td>0.60</td><td>0.76</td><td>6.58</td><td>0.31</td></tr><tr><td>A10 CC</td><td>3.25</td><td>1.13</td><td>1.22</td><td>5.49</td><td>0.77</td></tr><tr><td>AGN disturbed</td><td>3.29 ± 0.24</td><td>1.35 ± 0.09</td><td>1.38 ± 0.07</td><td>4.35 ± 0.16</td><td>0.39 ± 0.03]</td></tr><tr><td>P13 NCC</td><td>4.72</td><td>2.19</td><td>1.82</td><td>3.62</td><td>0.31</td></tr><tr><td>A10 dist.</td><td>3.20</td><td>1.08</td><td>1.41</td><td>5.49</td><td>0.38</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The median pressure profile of the MUSIC clusters (almost 1050) at redshift $z < 0.5$ for the AGN (left-hand panel), CSF (central panel), and NR (right-hand panel) simulation flavours. The grey lines represent the individual profiles for each cluster. The overall median profile is represented by black dots, with error bars computed by the MAD (equation 6). The yellow area shows the variation of the best fit for the gNFW between $1\sigma$ and $3\sigma$ intervals. The bottom panel of each plot shows the relative difference between the fit of the MUSIC profile and the corresponding profiles from A10 (green) and P13 (magenta), respectively.</figcaption></figure>

Piffaretti & Valdarnini 2008). The fit on the profile was performed in the radial range [$0.03–4 \times R_{500}$]. In this radial range, the observed profile is limited to $R_{500}$ and the region outside this radius was extrapolated according to the predictions from numerical simulations. They describe the resulting pressure profiles as universal, since it fits well both data from simulations and observations, with parameters listed in Table 2, fifth row.

The Planck Collaboration (Planck Collaboration V 2014a, hereafter P13) compared the median gNFW pressure profile of 62 nearby ($z < 0.5$) massive observed clusters with the A10 profile, finding that there is a very good agreement in the cluster intermediate radii between the two. However, within the core, i.e. $R < 0.15R_{500}$, the observed profile lies significantly below the A10 profile. The fit was done in the radial range [$0.02–3 \times R_{500}$] and the parameters are reported in the fourth row of Table 2.

To compare our simulated sample with observations, we compare the median pressure profiles (dots in Fig. 3) for all the MUSIC cluster sets at low redshifts ($z < 0.5$, a sample of almost 1050 clusters), with the A10 and P13 gNFW profiles, represented in green and violet, respectively. The three panels show the MUSIC

median profiles for the simulation flavours, AGN, CSF, and NR from left to right, with the MAD, see equation (6), as associated uncertainty. The yellow line and the shaded regions represent the MUSIC pressure profile fit, with the gNFW model parameters for AGN, CSF, and NR listed in the first three rows of Table 2. In the bottom plot of each panel, we present the relative difference between the MUSIC pressure profile fit, $f_M$, and the profile from A10 or P13

$$ \Delta = \frac{f_M - f_{A10/P13}}{f_M} \qquad (9) $$

The AGN flavour is the data set showing better agreement with both A10 and P13 profiles, especially for the A10’s universal profile, a very good approximation until $R_{500}$. In the case of the CSF flavour, the MUSIC profiles are steeper starting from $0.1R_{500}$, while in the NR case the MUSIC profiles are higher than the observed profiles within even a larger region (approximately $0.3R_{500}$). Since only the AGN set provides a reliable description of the observed profiles, we used this set to check the redshift dependence of the universal profile, extended to the high-redshift regime ($0.54 < z < 0.82$). As
<figure><img src="image_5.png" /><figcaption>Figure 4. Same as Fig. 3 but for the MUSIC clusters at high redshifts ($z > 0.5$, almost 800 clusters), for the AGN flavour only.</figcaption></figure>

we can see from Fig. 4, our data match well enough both the P13 profile and the A10 one with a relative difference of the order of 0.25. This difference is basically at the cluster core, inside $0.1R_{500}$, and in the outskirts, beyond $2R_{500}$. Contrary to the low $z$ case, now the P13 parameters are in better agreement with our simulations. The best-fitting parameters for the high-$z$ case are listed in Table 3.

We have also studied whether the dynamical state of the clusters could have an impact in the comparisons with observations. While we segregated extremely relaxed and disturbed clusters in MUSIC from the dynamical state parameters (see Section 4.3), it is not that straightforward with P13 observed data. The dynamical state classification of P13 is based simply on the presence or absence of a cool core (Planck Collaboration XI 2011b), and assumes that a cool core represents a relaxed system (see e.g. Hudson et al. 2010, although this assumption does not always hold, see Biffi et al. 2016). In A10, the cluster dynamical state classification is somewhat different: the classification was established in Pratt et al. (2009), where the subsamples were defined as cool core (10 systems), morphologically disturbed (12 systems), and neither cool core nor morphologically disturbed (11 systems). In Fig. 5, the median pressure profiles for relaxed and disturbed clusters are plotted both at low redshifts for the AGN flavour. The parameters of the fit are listed in Table 2, only for the AGN flavour because it is the one that better matches real observed clusters, as shown in Fig. 3. We can see that the disturbed profile from A10 and the NCC P13 profile match very well the MUSIC disturbed profile inside $R_{500}$. This means that the disturbed clusters from A10 and from the MUSIC simulation have similar features as the NCC clusters in P13. This does not happen for the relaxed MUSIC clusters, showing a shallower slope of the profile in the cluster core with respect to the CC clusters from A10 and P13. The reason for this behaviour could be that selecting dynamically relaxed clusters based on the indicators described in Section 4.3 does not allow us to discriminate CC clusters and clusters with a disturbed core. Therefore, we expect the inner slope of the MUSIC pressure profile estimated on relaxed clusters to be shallower than the one obtained by studying only CC clusters. Another possible explanation is that our AGN feedback model is too effective and expels more gas from the cluster core.

## 5.2 Electron density profile

The 3D electron density profiles, $n_e(r)$, of our simulated clusters are modelled using the analytical function proposed by Vikhlinin et al. (2006)

$$n_{\text{p}}n_{\text{e}}(r) = n_{0}^{2} \frac{(r/r_{\text{c}})^{-a}}{(1 + (r/r_{\text{c}})^{2})^{3b-a/2}} \frac{1}{(1 + (r/r_{\text{s}})^{c})^{e/c}} + \frac{n_{02}^{2}}{(1 + (r/r_{\text{c2}})^{2})^{3b_{2}}}, \quad (10)$$

which is a modification of the β-model (Cavaliere & Fusco-Femiano 1978) to represent the observed features of X-ray observations. It is based on two terms. The first term represents a cuspy profile near the cluster centre plus another power law to describe the steepening of the profile for $r > r_s$, the extra parameter c controls the width of the transition region. The second term is another β-model with a small central core to make the function flexible to fit the data in the central region of the clusters. All clusters profiles were fitted using the same parameter constraints as in Vikhlinin et al. (2006), that is fixing c = 3 and $e < 5$.

## 5.3 Temperature profile

The temperature profile is estimated as the mass-weighted average over an ensemble of gas particles within each spherical shell:

$$T = \frac{\sum_{i} (m_{i} \times T_{i})}{\sum_{i} m_{i}}, \quad (11)$$

where i run over each particle inside a spherical bin. We considered only particles with temperature $kT > 0.5$ keV to account only for the X-ray emitting gas. The analytical model describing the mass-weighted temperature was introduced also by Vikhlinin et al. (2006)

$$T(r) = T_0 \times \frac{x + \tau}{x + 1} \times \frac{(r/r_t)^{-a}}{(1 + (r/r_t)^b)^{c/b}} \\ x = (r/r_{\text{cool}})^{a_{\text{cool}}}. \quad (12)$$

The radial temperature profile has a broad peak at $r < 0.1R_{500}$ and decreases at larger radii, there is also a temperature decline towards the cluster centre, probably because of the presence of radiative cooling, represented by the central expression. Outside the central cooling region, the temperature profile is represented as a broken power law with a transition region, the last term, where $r_t$ is a scale radius. This model has eight free parameters, none of them was fixed.

The mass weighted temperature is the value that better relates to the mass of the cluster, in fact it reflects the gravitational potential of the system (Biffi et al. 2014). Nevertheless, there are various other ways of estimating the cluster temperature, for instance weighting the temperature by the X-ray emission of each gas particle, such as the spectroscopic-like temperature $T_{sl}$ (Mazzotta et al. 2004), in order to better explore the X-ray observable properties of simulated galaxy clusters and to compare against real observations (Rasia et al. 2014). Biffi et al. (2014) have computed the spectroscopic-like temperature for the MUSIC clusters by fitting the emission spectra from SPH particles in the energy band 0.5–10 keV. They found that this temperature is, on average, lower than the mass weighted temperature.

# 6 RESULTS FOR THE FULL MUSIC SAMPLE

In this section, we present the results for the hydrostatic masses defined by equations (1) and (2), for the complete MUSIC sample.
<figure><img src="image_6.png" /><figcaption>Figure 4. Same as Fig. 3 but for the MUSIC clusters at high redshifts ($z > 0.5$, almost 800 clusters), for the AGN flavour only.</figcaption></figure>

we can see from Fig. 4, our data match well enough both the P13 profile and the A10 one with a relative difference of the order of 0.25. This difference is basically at the cluster core, inside $0.1R_{500}$, and in the outskirts, beyond $2R_{500}$. Contrary to the low $z$ case, now the P13 parameters are in better agreement with our simulations. The best-fitting parameters for the high-$z$ case are listed in Table 3.

We have also studied whether the dynamical state of the clusters could have an impact in the comparisons with observations. While we segregated extremely relaxed and disturbed clusters in MUSIC from the dynamical state parameters (see Section 4.3), it is not that straightforward with P13 observed data. The dynamical state classification of P13 is based simply on the presence or absence of a cool core (Planck Collaboration XI 2011b), and assumes that a cool core represents a relaxed system (see e.g. Hudson et al. 2010, although this assumption does not always hold, see Biffi et al. 2016). In A10, the cluster dynamical state classification is somewhat different: the classification was established in Pratt et al. (2009), where the subsamples were defined as cool core (10 systems), morphologically disturbed (12 systems), and neither cool core nor morphologically disturbed (11 systems). In Fig. 5, the median pressure profiles for relaxed and disturbed clusters are plotted both at low redshifts for the AGN flavour. The parameters of the fit are listed in Table 2, only for the AGN flavour because it is the one that better matches real observed clusters, as shown in Fig. 3. We can see that the disturbed profile from A10 and the NCC P13 profile match very well the MUSIC disturbed profile inside $R_{500}$. This means that the disturbed clusters from A10 and from the MUSIC simulation have similar features as the NCC clusters in P13. This does not happen for the relaxed MUSIC clusters, showing a shallower slope of the profile in the cluster core with respect to the CC clusters from A10 and P13. The reason for this behaviour could be that selecting dynamically relaxed clusters based on the indicators described in Section 4.3 does not allow us to discriminate CC clusters and clusters with a disturbed core. Therefore, we expect the inner slope of the MUSIC pressure profile estimated on relaxed clusters to be shallower than the one obtained by studying only CC clusters. Another possible explanation is that our AGN feedback model is too effective and expels more gas from the cluster core.

## 5.2 Electron density profile

The 3D electron density profiles, $n_e(r)$, of our simulated clusters are modelled using the analytical function proposed by Vikhlinin et al. (2006)

$$n_{\text{p}}n_{\text{e}}(r) = n_{0}^{2} \frac{(r/r_{\text{c}})^{-a}}{(1 + (r/r_{\text{c}})^{2})^{3b-a/2}} \frac{1}{(1 + (r/r_{\text{s}})^{c})^{e/c}} + \frac{n_{02}^{2}}{(1 + (r/r_{\text{c2}})^{2})^{3b_{2}}}, \quad (10)$$

which is a modification of the β-model (Cavaliere & Fusco-Femiano 1978) to represent the observed features of X-ray observations. It is based on two terms. The first term represents a cuspy profile near the cluster centre plus another power law to describe the steepening of the profile for $r > r_s$, the extra parameter c controls the width of the transition region. The second term is another β-model with a small central core to make the function flexible to fit the data in the central region of the clusters. All clusters profiles were fitted using the same parameter constraints as in Vikhlinin et al. (2006), that is fixing c = 3 and e < 5.

## 5.3 Temperature profile

The temperature profile is estimated as the mass-weighted average over an ensemble of gas particles within each spherical shell:

$$T = \frac{\sum_{i} (m_{i} \times T_{i})}{\sum_{i} m_{i}}, \quad (11)$$

where i run over each particle inside a spherical bin. We considered only particles with temperature $kT > 0.5$ keV to account only for the X-ray emitting gas. The analytical model describing the mass-weighted temperature was introduced also by Vikhlinin et al. (2006)

$$T(r) = T_0 \times \frac{x + \tau}{x + 1} \times \frac{(r/r_t)^{-a}}{(1 + (r/r_t)^b)^{c/b}} \\ x = (r/r_{\text{cool}})^{a_{\text{cool}}}. \quad (12)$$

The radial temperature profile has a broad peak at $r < 0.1R_{500}$ and decreases at larger radii, there is also a temperature decline towards the cluster centre, probably because of the presence of radiative cooling, represented by the central expression. Outside the central cooling region, the temperature profile is represented as a broken power law with a transition region, the last term, where $r_t$ is a scale radius. This model has eight free parameters, none of them was fixed.

The mass weighted temperature is the value that better relates to the mass of the cluster, in fact it reflects the gravitational potential of the system (Biffi et al. 2014). Nevertheless, there are various other ways of estimating the cluster temperature, for instance weighting the temperature by the X-ray emission of each gas particle, such as the spectroscopic-like temperature $T_{sl}$ (Mazzotta et al. 2004), in order to better explore the X-ray observable properties of simulated galaxy clusters and to compare against real observations (Rasia et al. 2014). Biffi et al. (2014) have computed the spectroscopic-like temperature for the MUSIC clusters by fitting the emission spectra from SPH particles in the energy band 0.5–10 keV. They found that this temperature is, on average, lower than the mass weighted temperature.

# 6 RESULTS FOR THE FULL MUSIC SAMPLE

In this section, we present the results for the hydrostatic masses defined by equations (1) and (2), for the complete MUSIC sample.
<table><caption>**Table 3.** Same as Table 2 but for all the MUSIC clusters at high redshifts ($z > 0.5$) from AGN simulations (first row). The last two rows show the parameters when considering only relaxed and only disturbed clusters.</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>AGN</td><td>3.57 ± 0.13</td><td>1.62 ± 0.05</td><td>1.54 ± 0.05</td><td>4.18 ± 0.07</td><td>0.50 ± 0.01</td></tr><tr><td>AGN relaxed</td><td>4.52 ± 0.21</td><td>1.88 ± 0.05</td><td>1.61 ± 0.06</td><td>4.00 ± 0.07</td><td>0.53 ± 0.02</td></tr><tr><td>AGN disturbed</td><td>1.27 ± 0.17</td><td>1.44 ± 0.06</td><td>3.67 ± 0.74</td><td>3.93 ± 0.13</td><td>0.67 ± 0.05</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The best-fitting gNFW model for the median pressure profile of relaxed (in magenta) and disturbed (in blue) MUSIC clusters at low redshifts. The P13 profiles for CC and NCC clusters are also shown as dot-dashed lines (dark magenta and dark blue) and the A10 profiles as dashed lines. The bottom panel represents the relative difference, $\Delta$ (equation 9), for relaxed and disturbed clusters only.</figcaption></figure>

We will devote Section 7 and Appendix A to present the specific results for the NIKA2 twin sample. As we mentioned above, the HE mass estimates, and consequently, the mass bias, depend on different factors: the observable quantities used, the simulation flavour, the numerical method used to estimate the spatial gradients, and finally the redshift and the dynamical state of the clusters. In this section, we explore how the HE masses and biases depend on these factors.

## 6.1 Methods to estimate the HE masses

In this paper, we use two methods to estimate the HE masses accordingly to the ways of computing the derivatives in equations (1) and (2). In the first approach, the ICM radial profiles are estimated from the SPH gas particles within each spherical bin. Then the derivatives can be directly computed numerically from the binned profiles. This is a more direct estimation of the gradients but can also suffer from noise associated with the bin size and particle sampling. An alternative estimation is to first fit the numerical profiles by the analytical functions (described in Section 5) and take the derivative from the fitted function.

The first approach was applied in several studies (Ameglio et al. 2009; Lau et al. 2009; Battaglia et al. 2013; Sembolini et al. 2013; Nelson et al. 2014a; Biffi et al. 2016; Martizzi & Agrusa 2016; Shi et al. 2016; Le Brun et al. 2017; Cialone et al. 2018). In this work, each cluster is divided in spherical shells with 100 kpc width, which has been found to be the optimal binning for our purpose.

This procedure corresponds to smoothing the profile, still keeping its most important features when comparing with the original binning.

The second method is based on the analytical derivative of the fitting functions of the ICM radial profiles. Also in this case, several other works have made use of this method (Piffaretti & Valdarnini 2008; Meneghetti et al. 2010; Kay et al. 2012; Rasia et al. 2012; McCarthy et al. 2016; Gupta et al. 2017; Ruppin et al. 2019b; Ansarifard et al. 2020; Barnes et al. 2020; Pearce et al. 2020). The radial profiles of each cluster thermodynamic quantity are fitted with parametric models (see Sections 5.1, 5.2, and 5.3) using the PYTHON function *curve_fit*. The profile bootstrap errors have been estimated instead of the MADs, which mainly in the outskirts show too large values, often of the order of the median quantity. The fits are done in the radial range [0.1–1] × $R_{500}$. So, we neglect the core and the outskirts of the clusters, since in observed cores there is a large variation from cluster to cluster and it is still challenging to reach the external regions in X-ray. We fix a maximum value of 10 for the reduced $\chi^2$, hereafter $\hat{\chi}^2$, to include those fits which are still, visually, a good approximation to the profiles due to the slight variations and the small errors. In this way, we have almost 50 per cent of the total number of clusters which have reliable fits ($\hat{\chi}^2 \le 10$). The group of reliable fit clusters changes depending on the HE mass chosen (X-ray or SZ). In fact, having a cluster for a reliable fit for one of the thermodynamic quantities does not necessary mean that it has a good fit also for the other quantities.

## 6.2 Hydrostatic mass estimates

In the first row of Fig. 6, we show the HE masses computed at $R_{500}$, from SZ and X-ray observables estimated using the analytical fitting method, only for clusters that present reliable fits ($\hat{\chi}^2 < 10$) and for the AGN simulation flavour. They are plotted as a function of the cluster true mass $M_{500}$ at $z = 0$. The HE mass is proportional to the true mass. Clusters at other redshifts show a similar behaviour, as Le Brun et al. (2017) also find. We fit the linear relation between the HE mass and the true mass. The same analysis was done for all the simulation flavours and redshift bins. The values of the slope $a$, which is equivalent to the bias, $1 - b$, for all the simulation flavours and redshifts can be seen in Table 4. We do not find any clear dependence on redshift, with values of $1 - b$ that are around 0.8–0.9 for all the flavours, roughly in agreement also with the results using the numerical derivative method, not shown here. Taking into account only clusters with reliable fits leads to smaller slope values with respect to the numerical derivative method and to the situation in which we do not neglect any clusters with bad fits. In particular, in the last two methods usually the slope errors are larger than the only reliable fit method, and ultimately the results are compatible within $1\sigma$. However, fitting the profiles implies that we are not sensitive to extreme local ICM fluctuations, meaning that we are neglecting information about the clusters and that we could underestimate the
<table><caption>**Table 3.** Same as Table 2 but for all the MUSIC clusters at high redshifts ($z > 0.5$) from AGN simulations (first row). The last two rows show the parameters when considering only relaxed and only disturbed clusters.</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>AGN</td><td>3.57 ± 0.13</td><td>1.62 ± 0.05</td><td>1.54 ± 0.05</td><td>4.18 ± 0.07</td><td>0.50 ± 0.01</td></tr><tr><td>AGN relaxed</td><td>4.52 ± 0.21</td><td>1.88 ± 0.05</td><td>1.61 ± 0.06</td><td>4.00 ± 0.07</td><td>0.53 ± 0.02</td></tr><tr><td>AGN disturbed</td><td>1.27 ± 0.17</td><td>1.44 ± 0.06</td><td>3.67 ± 0.74</td><td>3.93 ± 0.13</td><td>0.67 ± 0.05</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The best-fitting gNFW model for the median pressure profile of relaxed (in magenta) and disturbed (in blue) MUSIC clusters at low redshifts. The P13 profiles for CC and NCC clusters are also shown as dot-dashed lines (dark magenta and dark blue) and the A10 profiles as dashed lines. The bottom panel represents the relative difference, $\Delta$ (equation 9), for relaxed and disturbed clusters only.</figcaption></figure>

We will devote Section 7 and Appendix A to present the specific results for the NIKA2 twin sample. As we mentioned above, the HE mass estimates, and consequently, the mass bias, depend on different factors: the observable quantities used, the simulation flavour, the numerical method used to estimate the spatial gradients, and finally the redshift and the dynamical state of the clusters. In this section, we explore how the HE masses and biases depend on these factors.

## 6.1 Methods to estimate the HE masses

In this paper, we use two methods to estimate the HE masses accordingly to the ways of computing the derivatives in equations (1) and (2). In the first approach, the ICM radial profiles are estimated from the SPH gas particles within each spherical bin. Then the derivatives can be directly computed numerically from the binned profiles. This is a more direct estimation of the gradients but can also suffer from noise associated with the bin size and particle sampling. An alternative estimation is to first fit the numerical profiles by the analytical functions (described in Section 5) and take the derivative from the fitted function.

The first approach was applied in several studies (Ameglio et al. 2009; Lau et al. 2009; Battaglia et al. 2013; Sembolini et al. 2013; Nelson et al. 2014a; Biffi et al. 2016; Martizzi & Agrusa 2016; Shi et al. 2016; Le Brun et al. 2017; Cialone et al. 2018). In this work, each cluster is divided in spherical shells with 100 kpc width, which has been found to be the optimal binning for our purpose.

This procedure corresponds to smoothing the profile, still keeping its most important features when comparing with the original binning.

The second method is based on the analytical derivative of the fitting functions of the ICM radial profiles. Also in this case, several other works have made use of this method (Piffaretti & Valdarnini 2008; Meneghetti et al. 2010; Kay et al. 2012; Rasia et al. 2012; McCarthy et al. 2016; Gupta et al. 2017; Ruppin et al. 2019b; Ansarifard et al. 2020; Barnes et al. 2020; Pearce et al. 2020). The radial profiles of each cluster thermodynamic quantity are fitted with parametric models (see Sections 5.1, 5.2, and 5.3) using the PYTHON function *curve_fit*. The profile bootstrap errors have been estimated instead of the MADs, which mainly in the outskirts show too large values, often of the order of the median quantity. The fits are done in the radial range [0.1–1] × $R_{500}$. So, we neglect the core and the outskirts of the clusters, since in observed cores there is a large variation from cluster to cluster and it is still challenging to reach the external regions in X-ray. We fix a maximum value of 10 for the reduced $\chi^2$, hereafter $\hat{\chi}^2$, to include those fits which are still, visually, a good approximation to the profiles due to the slight variations and the small errors. In this way, we have almost 50 per cent of the total number of clusters which have reliable fits ($\hat{\chi}^2 \le 10$). The group of reliable fit clusters changes depending on the HE mass chosen (X-ray or SZ). In fact, having a cluster for a reliable fit for one of the thermodynamic quantities does not necessary mean that it has a good fit also for the other quantities.

## 6.2 Hydrostatic mass estimates

In the first row of Fig. 6, we show the HE masses computed at $R_{500}$, from SZ and X-ray observables estimated using the analytical fitting method, only for clusters that present reliable fits ($\hat{\chi}^2 < 10$) and for the AGN simulation flavour. They are plotted as a function of the cluster true mass $M_{500}$ at $z = 0$. The HE mass is proportional to the true mass. Clusters at other redshifts show a similar behaviour, as Le Brun et al. (2017) also find. We fit the linear relation between the HE mass and the true mass. The same analysis was done for all the simulation flavours and redshift bins. The values of the slope $a$, which is equivalent to the bias, $1 - b$, for all the simulation flavours and redshifts can be seen in Table 4. We do not find any clear dependence on redshift, with values of $1 - b$ that are around 0.8–0.9 for all the flavours, roughly in agreement also with the results using the numerical derivative method, not shown here. Taking into account only clusters with reliable fits leads to smaller slope values with respect to the numerical derivative method and to the situation in which we do not neglect any clusters with bad fits. In particular, in the last two methods usually the slope errors are larger than the only reliable fit method, and ultimately the results are compatible within $1\sigma$. However, fitting the profiles implies that we are not sensitive to extreme local ICM fluctuations, meaning that we are neglecting information about the clusters and that we could underestimate the
<figure><img src="image_7.png" /></figure>

**Figure 6.** First row: The HE masses $M_{HE,SZ}$ and $M_{HE,X}$ as a function of $M_{500}$ for the AGN simulations at $z=0$ using the analytical fitting method. Only clusters with reliable fits are shown. A linear fit is performed and the value of the slopes are listed in the plots. The relaxed clusters are represented as magenta diamonds, the disturbed ones as blue squares, and the intermediate as black circles. Second row: A binned average of the ratio of $M_{HE}$ (the same of the first row) and $M_{500}$ is represented as a function of $M_{500}$. The magenta and blue lines are the averages of the relaxed and disturbed clusters, respectively, while the yellow one the average over all the clusters. The shaded regions show the $1\sigma$ intervals.

HE mass. On the other hand, gas substructures and fluctuations do influence the bias, possibly leading to an overestimation of it.

## 6.3 Hydrostatic mass biases

From the HE mass, the bias is estimated using equation (3). According to this definition, a mass bias value of 0 suggests an HE mass equal to the true mass, therefore the HE represents a good cluster mass approximation. This approximation does not account for contributions like, for instance, the non-thermal pressure (see Section 6.4).

At redshift 0, using the analytical fitting and the only reliable fits, we find, at $R_{500}$, $b_{SZ} = 0.23_{-0.09}^{+0.14}$ and $b_{X} = 0.14_{-0.13}^{+0.11}$ (represented in Fig. 1) for the AGN flavour. These results are given as median and 16th and 84th percentiles, since the bias distributions are not Gaussian, see Appendix B for more details. The biases at different redshifts and for different simulation flavour are listed in Table C1. We can see that the X-ray hydrostatic mass tends to give a better (closer) estimation of the real mass, with respect to the SZ one. This is mainly due to the way of computing the profiles of the gas, see Section 6.3.6, but the two are consistent within $1\sigma$.

### 6.3.1 Mass bias dependence on dynamical state

In order to study the link between the HE mass bias and the cluster dynamical state, we test the dependence on the relaxation parameter, $\chi_{DS}$ (see definition in 4.3). All of the clusters that have a $\chi_{DS} > 1$ are dynamically relaxed. The biases as a function of $\chi_{DS}$, at redshift 0, are shown in Fig. 7 for the analytical fitting method, only for reliable

fits. The AGN flavour is shown here, but the situation is similar for CSF and NR. The relaxed clusters present a lower scatter, compared to the other cases and also a smaller overall HE bias, as it is shown in this plot and in Table C1. For the SZ and X-ray bias, the Pearson correlation coefficient is, respectively, -0.4 and -0.3, resulting in a moderate correlation and the slope of the linear relation between the bias and $\chi_{DS}$ is smaller and different from zero at $3\sigma$ level. The latter observation is true also for the numerical derivative method (which employs all the clusters of the sample), yet the correlation is weaker, around -0.2 for both biases.

### 6.3.2 Mass bias dependence on radial profile

The numerical derivative method is powerful to explore the HE mass bias even in the large radial range. The median radial profiles and the MADs of $b_{SZ}$ and $b_{X}$ at $z=0$ are shown in Fig. 8, from 0.2 to 2 $R_{500}$. In the left-hand panel of each row, there is the median profile over all the clusters, in the middle only relaxed, and in the right only disturbed clusters. The bias profiles from all and only relaxed classes are very regular and, as expected, increase at large radii, as found also in other works (e.g. Meneghetti et al. 2010; Cialone et al. 2018; Ansarifard et al. 2020). This could be most likely the presence of higher non-thermal processes in the cluster outskirts, see Section 6.4. On the contrary, the trend of the disturbed clusters has a lot of variations and a large scatter, with a dip between $R_{500}$ and $1.5R_{500}$, which seems not due to the non-thermal pressure, as a matter of fact it is still present after applying the non-thermal correction, see Section 6.4. We will study in deep this behaviour in a future work, but we attempt here a possible explanation. It could be due to merger shocks that could boost either the ICM temperature (and therefore $M_{HE}$) and consequently the gas pressure.

### 6.3.3 Bias dependence on baryon models

Taking advantage of exploring the same objects with different flavours, we are able to compare HE mass bias values, without worrying about the simulation features, like resolution, integrating box size, cluster mass range, and number of objects, focusing only on the differences due to the physics included in the simulation. At redshift 0 and at $R_{500}$, using the analytical fitting and the only reliable fits, we find $b_{SZ} = 0.23_{-0.09}^{+0.14}$ and $b_{X} = 0.14_{-0.13}^{+0.11}$ for the AGN flavour, $b_{SZ} = 0.26_{-0.10}^{+0.12}$ and $b_{X} = 0.14_{-0.11}^{+0.16}$ for the CSF flavour, and $b_{SZ} = 0.27_{-0.11}^{+0.15}$ and $b_{X} = 0.15_{-0.13}^{+0.14}$ for NR. As we can see, the simulation flavour does not influence much the bias value, as it is clear also in Table C1. In Fig. 8, the radial profiles of both biases are represented for each flavour. Also in this case, the all clusters median profiles are not affected by the baryonic physics. From the middle and right-hand panels, we see that this trend is common also to the relaxed clusters while for the disturbed objects there is a tendency for the AGN runs to have a smaller bias, although the profiles are all consistent within the scatter.

### 6.3.4 The bias dependence on redshift

The bias dependence on the redshift is presented in the first row of Fig. 9. The results are from the analytical fitting method with only the reliable fits and the AGN flavour (in the case of the other flavours the situation is similar). With the current errors, we do not detect any dependence with the redshift, for either the bias (as in Le Brun et al. 2017; Henson et al. 2017; Salvati et al. 2018; Koukoufilippas et al. 2020) or its uncertainties. In all cases, the disturbed clusters
<figure><img src="image_8.png" /></figure>

**Figure 6.** First row: The HE masses $M_{HE,SZ}$ and $M_{HE,X}$ as a function of $M_{500}$ for the AGN simulations at $z=0$ using the analytical fitting method. Only clusters with reliable fits are shown. A linear fit is performed and the value of the slopes are listed in the plots. The relaxed clusters are represented as magenta diamonds, the disturbed ones as blue squares, and the intermediate as black circles. Second row: A binned average of the ratio of $M_{HE}$ (the same of the first row) and $M_{500}$ is represented as a function of $M_{500}$. The magenta and blue lines are the averages of the relaxed and disturbed clusters, respectively, while the yellow one the average over all the clusters. The shaded regions show the $1\sigma$ intervals.

HE mass. On the other hand, gas substructures and fluctuations do influence the bias, possibly leading to an overestimation of it.

## 6.3 Hydrostatic mass biases

From the HE mass, the bias is estimated using equation (3). According to this definition, a mass bias value of 0 suggests an HE mass equal to the true mass, therefore the HE represents a good cluster mass approximation. This approximation does not account for contributions like, for instance, the non-thermal pressure (see Section 6.4).

At redshift 0, using the analytical fitting and the only reliable fits, we find, at $R_{500}$, $b_{SZ} = 0.23_{-0.09}^{+0.14}$ and $b_X = 0.14_{-0.13}^{+0.11}$ (represented in Fig. 1) for the AGN flavour. These results are given as median and 16th and 84th percentiles, since the bias distributions are not Gaussian, see Appendix B for more details. The biases at different redshifts and for different simulation flavour are listed in Table C1. We can see that the X-ray hydrostatic mass tends to give a better (closer) estimation of the real mass, with respect to the SZ one. This is mainly due to the way of computing the profiles of the gas, see Section 6.3.6, but the two are consistent within $1\sigma$.

### 6.3.1 Mass bias dependence on dynamical state

In order to study the link between the HE mass bias and the cluster dynamical state, we test the dependence on the relaxation parameter, $\chi_{DS}$ (see definition in 4.3). All of the clusters that have a $\chi_{DS} > 1$ are dynamically relaxed. The biases as a function of $\chi_{DS}$, at redshift 0, are shown in Fig. 7 for the analytical fitting method, only for reliable

fits. The AGN flavour is shown here, but the situation is similar for CSF and NR. The relaxed clusters present a lower scatter, compared to the other cases and also a smaller overall HE bias, as it is shown in this plot and in Table C1. For the SZ and X-ray bias, the Pearson correlation coefficient is, respectively, -0.4 and -0.3, resulting in a moderate correlation and the slope of the linear relation between the bias and $\chi_{DS}$ is smaller and different from zero at $3\sigma$ level. The latter observation is true also for the numerical derivative method (which employs all the clusters of the sample), yet the correlation is weaker, around -0.2 for both biases.

### 6.3.2 Mass bias dependence on radial profile

The numerical derivative method is powerful to explore the HE mass bias even in the large radial range. The median radial profiles and the MADs of $b_{SZ}$ and $b_X$ at $z=0$ are shown in Fig. 8, from 0.2 to 2 $R_{500}$. In the left-hand panel of each row, there is the median profile over all the clusters, in the middle only relaxed, and in the right only disturbed clusters. The bias profiles from all and only relaxed classes are very regular and, as expected, increase at large radii, as found also in other works (e.g. Meneghetti et al. 2010; Cialone et al. 2018; Ansarifard et al. 2020). This could be most likely the presence of higher non-thermal processes in the cluster outskirts, see Section 6.4. On the contrary, the trend of the disturbed clusters has a lot of variations and a large scatter, with a dip between $R_{500}$ and $1.5R_{500}$, which seems not due to the non-thermal pressure, as a matter of fact it is still present after applying the non-thermal correction, see Section 6.4. We will study in deep this behaviour in a future work, but we attempt here a possible explanation. It could be due to merger shocks that could boost either the ICM temperature (and therefore $M_{HE}$) and consequently the gas pressure.

### 6.3.3 Bias dependence on baryon models

Taking advantage of exploring the same objects with different flavours, we are able to compare HE mass bias values, without worrying about the simulation features, like resolution, integrating box size, cluster mass range, and number of objects, focusing only on the differences due to the physics included in the simulation. At redshift 0 and at $R_{500}$, using the analytical fitting and the only reliable fits, we find $b_{SZ} = 0.23_{-0.09}^{+0.14}$ and $b_X = 0.14_{-0.13}^{+0.11}$ for the AGN flavour, $b_{SZ} = 0.26_{-0.10}^{+0.12}$ and $b_X = 0.14_{-0.11}^{+0.16}$ for the CSF flavour, and $b_{SZ} = 0.27_{-0.11}^{+0.15}$ and $b_X = 0.15_{-0.13}^{+0.14}$ for NR. As we can see, the simulation flavour does not influence much the bias value, as it is clear also in Table C1. In Fig. 8, the radial profiles of both biases are represented for each flavour. Also in this case, the all clusters median profiles are not affected by the baryonic physics. From the middle and right-hand panels, we see that this trend is common also to the relaxed clusters while for the disturbed objects there is a tendency for the AGN runs to have a smaller bias, although the profiles are all consistent within the scatter.

### 6.3.4 The bias dependence on redshift

The bias dependence on the redshift is presented in the first row of Fig. 9. The results are from the analytical fitting method with only the reliable fits and the AGN flavour (in the case of the other flavours the situation is similar). With the current errors, we do not detect any dependence with the redshift, for either the bias (as in Le Brun et al. 2017; Henson et al. 2017; Salvati et al. 2018; Koukoufilippas et al. 2020) or its uncertainties. In all cases, the disturbed clusters
<table><caption>**Table 4.** Slopes and $1\sigma$ errors from the fit $M_{HE} = aM_{500}$ in Fig. 6 at all the redshifts and for all the simulation flavour with the analytical fitting method using only the clusters with reliable fits.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Slope a<sub>SZ</sub></th><th colspan="3">Slope a<sub>X</sub></th></tr><tr><th>AGN</th><th>CSF</th><th>NR</th><th>AGN</th><th>CSF</th><th>NR</th></tr></thead><tbody><tr><td>0.0</td><td>0.808 ± 0.014</td><td>0.790 ± 0.013</td><td>0.763 ± 0.016</td><td>0.916 ± 0.016</td><td>0.909 ± 0.015</td><td>0.890 ± 0.018</td></tr><tr><td>0.11</td><td>0.819 ± 0.014</td><td>0.800 ± 0.014</td><td>0.767 ± 0.016</td><td>0.914 ± 0.012</td><td>0.922 ± 0.015</td><td>0.897 ± 0.014</td></tr><tr><td>0.33</td><td>0.788 ± 0.015</td><td>0.787 ± 0.013</td><td>0.731 ± 0.010</td><td>0.876 ± 0.015</td><td>0.860 ± 0.013</td><td>0.871 ± 0.018</td></tr><tr><td>0.43</td><td>0.783 ± 0.016</td><td>0.766 ± 0.010</td><td>0.716 ± 0.011</td><td>0.907 ± 0.017</td><td>0.886 ± 0.013</td><td>0.889 ± 0.018</td></tr><tr><td>0.54</td><td>0.758 ± 0.014</td><td>0.739 ± 0.012</td><td>0.758 ± 0.015</td><td>0.880 ± 0.015</td><td>0.886 ± 0.011</td><td>0.855 ± 0.014</td></tr><tr><td>0.67</td><td>0.778 ± 0.013</td><td>0.752 ± 0.009</td><td>0.764 ± 0.016</td><td>0.898 ± 0.017</td><td>0.923 ± 0.013</td><td>0.983 ± 0.026</td></tr><tr><td>0.82</td><td>0.763 ± 0.014</td><td>0.789 ± 0.014</td><td>0.712 ± 0.010</td><td>0.902 ± 0.018</td><td>0.895 ± 0.015</td><td>0.902 ± 0.021</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** The dependence of the SZ bias (left) and the X-ray bias (right) on the relaxation parameter $\chi_{DS}$. The vertical dot-dashed line represents the threshold dividing the dynamically relaxed cluster (right) from disturbed ones (left). In the legend, the best-fitting (blue solid line) values and their errors are shown. Following the classification of the 3D morphological parameters, Section 4.3, the relaxed clusters are represented as magenta diamonds, the blue squares are the disturbed ones, and the intermediate corresponds to black circles.</figcaption></figure>

have the largest errors, as already pointed out in Section 6.3.1. We see the same behaviour for the numerical derivative method, with the only difference that the disturbed cluster percentiles are almost twice larger than the percentiles shown in the analytical fitting method. In this figure, the bootstrap error is represented too, for the three cases.

### 6.3.5 The bias mass dependence

A binned average of the quantity $1 - b = M_{HE}/M_{500}$ is represented in the second row of Fig. 6, as a function of $M_{500}$. We do not find any mass dependence inside the $1\sigma$ regions, in agreement with Ansarifard et al. (2020).

Several authors detect a dependence when using the spectroscopic-like temperature, see e.g. Rasia et al. (2006), Piffaretti & Valdarnini (2008), Le Brun et al. (2017), Henson et al. (2017), Pearce et al. (2020), Barnes et al. (2020). This result is mainly due to the fitting of the gas multitemperature spectrum with a single-temperature model. The plasma in the outskirts is characterized by a higher degree of substructures not yet thermalized and therefore colder and denser. This reduces the temperature evaluated at $T_{500}$. Indeed, the mass bias estimated using the spectroscopic-like temperature, is generally higher and with a larger scatter with respect to the mass-weighted temperature case. We have also confirmed this behaviour in our MUSIC HE bias estimates (at redshift 0) when using the spectroscopic-like temperatures, mostly because of the larger scatter of their profiles. But, again, no dependence of the HE bias on the cluster mass was found when the spectroscopic-like temperature profiles are used.

### 6.3.6 Correlation between b<sub>X</sub> and b<sub>SZ</sub>

We find that the two bias factors are strongly correlated, with a Pearson coefficient of 0.8. To estimate the scatter between the two biases, we can use the mean absolute difference

$$MD = \sqrt{\frac{\sum_i (b_{X,i} - b_{SZ,i})^2}{N(N-1)}} \quad (13)$$

with $N$ being the total number of clusters. At all redshifts $MD$ is of the order of $10^{-1}$.

In the simulations, the pressure of each fluid particle is computed assuming the ideal gas equation of state, used also to derive $M_{HE,X}$, equation (2), from $M_{HE,SZ}$, equation (1),

$$P_i = n_{e,i} T_i, \quad (14)$$

where $n_{e,i}$ and $T_i$ are the electron density and the temperature of each particle of fluid. Then, the median over the spherical shell is performed, obtaining the median pressure profile. This profile, together with the median profile of the electron density and the mass-weighted temperature profile, are used in equations (1) and (2). For this reason, we expect the two HE masses to be different.

To verify that the scatter between the two biases is really due to the previous consideration, we need to quantify it, for instance using

$$\frac{M_{HE,SZ}}{M_{HE,X}} = \frac{1 - b_{SZ}}{1 - b_{X}}, \quad (15)$$

which, replacing the masses with the equations (1) and (2), becomes

$$\frac{1 - b_{SZ}}{1 - b_{X}} = \frac{P}{n_{e}T} \frac{d \ln P / d \ln r}{d \ln(n_{e}T) / d \ln r} \sim \frac{P}{n_{e}T}, \quad (16)$$

where the derivatives fraction is always of the order of 1, no matter the flavour or the redshift.

We find that the profile of $P/n_{e}T$ decreases very slowly from a value of almost 1 in the core. In the case of AGN, $z = 0$, the deviation from 1 of this ratio is $1 - P/n_{e}T = 0.12$, but in general it is always of the order of $10^{-1}$. This is of the same order of the scatter $MD$, so we can conclude that the difference between the X-ray and SZ bias is mainly due to the use of the median and mass-weighted profiles.

## 6.4 The non-thermal pressure contribution

The HE approximation does not take into account non-thermal motions of the ICM, which could have a significant contribution to the pressure support of the gas within the cluster gravitational potential. The non-thermal pressure contribution comes mostly from turbulence or bulk motions of the ICM, neglecting other sources such as magnetic fields or cosmic rays pressure. Recent studies have shown that it can contribute as much as 30 per cent or more to the overall gas pressure at $R_{500}$ (Lau et al. 2009; Nelson et al. 2014a;
<table><caption>**Table 4.** Slopes and $1\sigma$ errors from the fit $M_{HE} = aM_{500}$ in Fig. 6 at all the redshifts and for all the simulation flavour with the analytical fitting method using only the clusters with reliable fits.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Slope a<sub>SZ</sub></th><th colspan="3">Slope a<sub>X</sub></th></tr><tr><th>AGN</th><th>CSF</th><th>NR</th><th>AGN</th><th>CSF</th><th>NR</th></tr></thead><tbody><tr><td>0.0</td><td>0.808 ± 0.014</td><td>0.790 ± 0.013</td><td>0.763 ± 0.016</td><td>0.916 ± 0.016</td><td>0.909 ± 0.015</td><td>0.890 ± 0.018</td></tr><tr><td>0.11</td><td>0.819 ± 0.014</td><td>0.800 ± 0.014</td><td>0.767 ± 0.016</td><td>0.914 ± 0.012</td><td>0.922 ± 0.015</td><td>0.897 ± 0.014</td></tr><tr><td>0.33</td><td>0.788 ± 0.015</td><td>0.787 ± 0.013</td><td>0.731 ± 0.010</td><td>0.876 ± 0.015</td><td>0.860 ± 0.013</td><td>0.871 ± 0.018</td></tr><tr><td>0.43</td><td>0.783 ± 0.016</td><td>0.766 ± 0.010</td><td>0.716 ± 0.011</td><td>0.907 ± 0.017</td><td>0.886 ± 0.013</td><td>0.889 ± 0.018</td></tr><tr><td>0.54</td><td>0.758 ± 0.014</td><td>0.739 ± 0.012</td><td>0.758 ± 0.015</td><td>0.880 ± 0.015</td><td>0.886 ± 0.011</td><td>0.855 ± 0.014</td></tr><tr><td>0.67</td><td>0.778 ± 0.013</td><td>0.752 ± 0.009</td><td>0.764 ± 0.016</td><td>0.898 ± 0.017</td><td>0.923 ± 0.013</td><td>0.983 ± 0.026</td></tr><tr><td>0.82</td><td>0.763 ± 0.014</td><td>0.789 ± 0.014</td><td>0.712 ± 0.010</td><td>0.902 ± 0.018</td><td>0.895 ± 0.015</td><td>0.902 ± 0.021</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The dependence of the SZ bias (left) and the X-ray bias (right) on the relaxation parameter $\chi_{DS}$. The vertical dot-dashed line represents the threshold dividing the dynamically relaxed cluster (right) from disturbed ones (left). In the legend, the best-fitting (blue solid line) values and their errors are shown. Following the classification of the 3D morphological parameters, Section 4.3, the relaxed clusters are represented as magenta diamonds, the blue squares are the disturbed ones, and the intermediate corresponds to black circles.</figcaption></figure>

have the largest errors, as already pointed out in Section 6.3.1. We see the same behaviour for the numerical derivative method, with the only difference that the disturbed cluster percentiles are almost twice larger than the percentiles shown in the analytical fitting method. In this figure, the bootstrap error is represented too, for the three cases.

### 6.3.5 The bias mass dependence

A binned average of the quantity $1 - b = M_{HE}/M_{500}$ is represented in the second row of Fig. 6, as a function of $M_{500}$. We do not find any mass dependence inside the $1\sigma$ regions, in agreement with Ansarifard et al. (2020).

Several authors detect a dependence when using the spectroscopic-like temperature, see e.g. Rasia et al. (2006), Piffaretti & Valdarnini (2008), Le Brun et al. (2017), Henson et al. (2017), Pearce et al. (2020), Barnes et al. (2020). This result is mainly due to the fitting of the gas multitemperature spectrum with a single-temperature model. The plasma in the outskirts is characterized by a higher degree of substructures not yet thermalized and therefore colder and denser. This reduces the temperature evaluated at $T_{500}$. Indeed, the mass bias estimated using the spectroscopic-like temperature, is generally higher and with a larger scatter with respect to the mass-weighted temperature case. We have also confirmed this behaviour in our MUSIC HE bias estimates (at redshift 0) when using the spectroscopic-like temperatures, mostly because of the larger scatter of their profiles. But, again, no dependence of the HE bias on the cluster mass was found when the spectroscopic-like temperature profiles are used.

### 6.3.6 Correlation between b<sub>X</sub> and b<sub>SZ</sub>

We find that the two bias factors are strongly correlated, with a Pearson coefficient of 0.8. To estimate the scatter between the two biases, we can use the mean absolute difference

$$MD = \sqrt{\frac{\sum_i (b_{X,i} - b_{SZ,i})^2}{N(N-1)}} \quad (13)$$

with $N$ being the total number of clusters. At all redshifts MD is of the order of $10^{-1}$.

In the simulations, the pressure of each fluid particle is computed assuming the ideal gas equation of state, used also to derive $M_{HE,X}$, equation (2), from $M_{HE,SZ}$, equation (1),

$$P_i = n_{e,i} T_i, \quad (14)$$

where $n_{e,i}$ and $T_i$ are the electron density and the temperature of each particle of fluid. Then, the median over the spherical shell is performed, obtaining the median pressure profile. This profile, together with the median profile of the electron density and the mass-weighted temperature profile, are used in equations (1) and (2). For this reason, we expect the two HE masses to be different.

To verify that the scatter between the two biases is really due to the previous consideration, we need to quantify it, for instance using

$$\frac{M_{HE,SZ}}{M_{HE,X}} = \frac{1 - b_{SZ}}{1 - b_{X}}, \quad (15)$$

which, replacing the masses with the equations (1) and (2), becomes

$$\frac{1 - b_{SZ}}{1 - b_{X}} = \frac{P}{n_{e}T} \frac{d \ln P / d \ln r}{d \ln(n_{e}T) / d \ln r} \sim \frac{P}{n_{e}T}, \quad (16)$$

where the derivatives fraction is always of the order of 1, no matter the flavour or the redshift.

We find that the profile of $P/n_{e}T$ decreases very slowly from a value of almost 1 in the core. In the case of AGN, $z = 0$, the deviation from 1 of this ratio is $1 - P/n_{e}T = 0.12$, but in general it is always of the order of $10^{-1}$. This is of the same order of the scatter MD, so we can conclude that the difference between the X-ray and SZ bias is mainly due to the use of the median and mass-weighted profiles.

## 6.4 The non-thermal pressure contribution

The HE approximation does not take into account non-thermal motions of the ICM, which could have a significant contribution to the pressure support of the gas within the cluster gravitational potential. The non-thermal pressure contribution comes mostly from turbulence or bulk motions of the ICM, neglecting other sources such as magnetic fields or cosmic rays pressure. Recent studies have shown that it can contribute as much as 30 per cent or more to the overall gas pressure at $R_{500}$ (Lau et al. 2009; Nelson et al. 2014a;
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The median SZ bias radial profile (upper row) and the corresponding one for the X-ray bias (bottom row) for $z = 0$ and the three simulated flavours. In each panel, the median profile is plotted with a blue line for AGN, the CSF profile with a green line, and the NR with a red line. The shaded regions represent the 1 MAD interval for each case. In the left-hand panels, the bias is estimated over all the MUSIC clusters, in the middle, over only the relaxed ones, and in the right-hand panel only over the disturbed clusters.</figcaption></figure>

Biffi et al. 2016; Martizzi & Agrusa 2016; Shi et al. 2016; Vazza et al. 2018; Pearce et al. 2020). This has been identified as the main origin of the cluster mass bias.

The non-thermal pressure component is modelled as

$$P_{\text{nth}} = \alpha(r) P_{\text{tot}} \qquad (17)$$

assuming $P_{\text{tot}} = P_{\text{th}} + P_{\text{nth}}$. From simulations, the non-thermal pressure (Nelson, Lau & Nagai 2014b) can be estimated as

$$P_{\text{nth}} = \rho \sigma^2, \qquad (18)$$

where $\sigma^2 = \sum_j \sigma_j^2$, with $j = (x, y, z)$, is the square of the three-dimensional velocity dispersion of gas particles and

$$\sigma_j^2 = \frac{\sum_i (m_i v_{i,j} - \bar{v}_j)^2}{\sum_i m_i}, \qquad (19)$$

is the velocity dispersion in each spatial direction computed from all the gas particles within each spherical shell around the cluster centre.

Then, introducing $P_{\text{tot}}$ in the formula for the HE, we can derive the corrected hydrostatic masses (Pearce et al. 2020)

$$M_{\text{HE,SZ,corr}} = \frac{1}{1-\alpha} \left[ M_{\text{HE,SZ}} - \frac{\alpha}{1-\alpha} \frac{r P_{\text{th}}}{G \mu m_{\rho} n_{e}} \frac{d \ln \alpha}{d \ln r} \right] \qquad (20)$$

$$M_{\text{HE,X,corr}} = \frac{1}{1-\alpha} \left[ M_{\text{X,SZ}} - \frac{\alpha}{1-\alpha} \frac{k_{\text{B}} T r}{G \mu m_{\rho}} \frac{d \ln \alpha}{d \ln r} \right]. \qquad (21)$$

We note that the temperature profile used in the formula above corresponds to the mass-weighted temperature definition. For the non-thermal correction applied on HE mass using the spectroscopic-like temperature, see Rasia et al. (2006).

Therefore, in order to correct the HE masses we need to know the value of the $\alpha(r)$ function. We use the fitting formula originally proposed by Nelson et al. (2014b) for $R_{200}$, adapted to the smaller aperture radius $R_{500}$,

$$\alpha(r) = 1 - A \left[ 1 + \exp \left( - \left( \frac{r/R_{500}}{B} \right)^C \right) \right] \qquad (22)$$

where A, B, and C are three free parameters (Pearce et al. 2020).

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The redshift evolution of the biases, $b_{\text{SZ}}$ (left-hand panels) and $b_{\text{X}}$ (right-hand panels), from the AGN simulations. The reliable fits have been used to compute them. The median values of the bias for all clusters, relaxed, and disturbed are represented with yellow triangles, magenta diamonds, and blue squares, respectively. The error bars correspond to bootstrap errors. The shaded regions represent the 16th and 84th percentiles. In the bottom panels, the corrected bias factors $b_{\text{corr}}$ (see Section 6.4) are also plotted.</figcaption></figure>

The contribution of non-thermal pressure becomes significant at large radii (e.g. Biffi et al. 2016; Pearce et al. 2020). We confirm a similar behaviour as shown in Fig. 10, where the ratio $\alpha$ is plotted. We can see that, at $R_{500}$, it goes from ~15 per cent to ~40 per cent, depending on the morphological state of the clusters. The most disturbed clusters have an higher contribution to $\alpha$, even if the two classes of objects (relaxed and disturbed) are consistent within their errors. Applying the non-thermal correction to the mass should result in minimizing the biases. In the bottom panels of Fig. 9, indeed, the $b_{\text{X}}$ and $b_{\text{SZ}}$ are close to 0, but the scatter (represented by the percentiles in Table C1) is larger than the non-corrected case due to the large scatter in the $\alpha$ profiles, in agreement with the results from Pearce et al. (2020). In addition, we also show that there is not a significant
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The median SZ bias radial profile (upper row) and the corresponding one for the X-ray bias (bottom row) for $z = 0$ and the three simulated flavours. In each panel, the median profile is plotted with a blue line for AGN, the CSF profile with a green line, and the NR with a red line. The shaded regions represent the 1 MAD interval for each case. In the left-hand panels, the bias is estimated over all the MUSIC clusters, in the middle, over only the relaxed ones, and in the right-hand panel only over the disturbed clusters.</figcaption></figure>

Biffi et al. 2016; Martizzi & Agrusa 2016; Shi et al. 2016; Vazza et al. 2018; Pearce et al. 2020). This has been identified as the main origin of the cluster mass bias.

The non-thermal pressure component is modelled as

$$P_{\text{nth}} = \alpha(r) P_{\text{tot}} \qquad (17)$$

assuming $P_{\text{tot}} = P_{\text{th}} + P_{\text{nth}}$. From simulations, the non-thermal pressure (Nelson, Lau & Nagai 2014b) can be estimated as

$$P_{\text{nth}} = \rho \sigma^2, \qquad (18)$$

where $\sigma^2 = \sum_j \sigma_j^2$, with $j = (x, y, z)$, is the square of the three-dimensional velocity dispersion of gas particles and

$$\sigma_j^2 = \frac{\sum_i (m_i v_{i,j} - \bar{v}_j)^2}{\sum_i m_i}, \qquad (19)$$

is the velocity dispersion in each spatial direction computed from all the gas particles within each spherical shell around the cluster centre.

Then, introducing $P_{\text{tot}}$ in the formula for the HE, we can derive the corrected hydrostatic masses (Pearce et al. 2020)

$$M_{\text{HE,SZ,corr}} = \frac{1}{1-\alpha} \left[ M_{\text{HE,SZ}} - \frac{\alpha}{1-\alpha} \frac{r P_{\text{th}}}{G \mu m_{\rho} n_{e}} \frac{d \ln \alpha}{d \ln r} \right] \qquad (20)$$

$$M_{\text{HE,X,corr}} = \frac{1}{1-\alpha} \left[ M_{\text{X,SZ}} - \frac{\alpha}{1-\alpha} \frac{k_{\text{B}} T r}{G \mu m_{\rho}} \frac{d \ln \alpha}{d \ln r} \right]. \qquad (21)$$

We note that the temperature profile used in the formula above corresponds to the mass-weighted temperature definition. For the non-thermal correction applied on HE mass using the spectroscopic-like temperature, see Rasia et al. (2006).

Therefore, in order to correct the HE masses we need to know the value of the $\alpha(r)$ function. We use the fitting formula originally proposed by Nelson et al. (2014b) for $R_{200}$, adapted to the smaller aperture radius $R_{500}$,

$$\alpha(r) = 1 - A \left[ 1 + \exp \left( - \left( \frac{r/R_{500}}{B} \right)^C \right) \right] \qquad (22)$$

where A, B, and C are three free parameters (Pearce et al. 2020).

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The redshift evolution of the biases, $b_{\text{SZ}}$ (left-hand panels) and $b_{\text{X}}$ (right-hand panels), from the AGN simulations. The reliable fits have been used to compute them. The median values of the bias for all clusters, relaxed, and disturbed are represented with yellow triangles, magenta diamonds, and blue squares, respectively. The error bars correspond to bootstrap errors. The shaded regions represent the 16th and 84th percentiles. In the bottom panels, the corrected bias factors $b_{\text{corr}}$ (see Section 6.4) are also plotted.</figcaption></figure>

The contribution of non-thermal pressure becomes significant at large radii (e.g. Biffi et al. 2016; Pearce et al. 2020). We confirm a similar behaviour as shown in Fig. 10, where the ratio $\alpha$ is plotted. We can see that, at $R_{500}$, it goes from ~15 per cent to ~40 per cent, depending on the morphological state of the clusters. The most disturbed clusters have an higher contribution to $\alpha$, even if the two classes of objects (relaxed and disturbed) are consistent within their errors. Applying the non-thermal correction to the mass should result in minimizing the biases. In the bottom panels of Fig. 9, indeed, the $b_{\text{X}}$ and $b_{\text{SZ}}$ are close to 0, but the scatter (represented by the percentiles in Table C1) is larger than the non-corrected case due to the large scatter in the $\alpha$ profiles, in agreement with the results from Pearce et al. (2020). In addition, we also show that there is not a significant
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The mean radial profile of $\alpha = P_{\text{nth}}/P_{\text{tot}}$ for MUSIC clusters in the AGN simulation at $z=0$. The profiles for all the clusters, the only relaxed, and the only disturbed are shown as yellow, magenta, and blue solid lines, respectively. The shaded regions represent the $1\sigma$ interval.</figcaption></figure>

dependence on the redshift range analysed in this work. Moreover, the correction on $b_X$, which already gave an estimation of the mass nearer to 0 before the correction, is more effective than on the $b_{SZ}$. This leads to zero X-ray bias for the relaxed and all cluster cases, and to an overcorrection for the disturbed class. Therefore, the dynamical state has an impact on the strength of the correction, but not on its radial profile, as stated before and shown in Fig. 10.

The median radial profiles of the two corrected biases $b_{SZ,corr}$ and $b_{X,corr}$ are presented in Fig. 11 for redshift 0 and all flavours. The increasing radial dependence of the X-ray and SZ bias was cancelled when including the non-thermal correction, as expected. However, there is still a large variation in the bias profiles around $R_{500}$ for the disturbed clusters, which was already present before the correction. This could be explained by a more intriguing effect, which does not depend on the non-thermal pressure, or the simple correction formula does not account for the extreme non-thermal pressure from the disturbed clusters. We also note that deviations from spherical symmetry as well as the presence of substructures impacting the multiphase structure of the gas can play a significant role in the outskirts. As we pointed out before, the scatter (represented by the MADs in this plot) is larger than the non-corrected bias. Moreover, the SZ bias has larger scatter than the X-ray bias, especially for the CSF and NR flavours. The NR biases have less regular profiles with respect to the other flavours, but there is not a substantial difference between them. All the different flavour profiles are in agreement.

# 7 RESULTS FROM THE NIKA2 TWIN SAMPLE

The NIKA2 twin sample (Ruppin et al. 2019b) is composed of 32 MUSIC clusters, see Section 4.2. To build the sample, we chose clusters at $z = 0.82$ that are not the progenitors of the clusters at $z = 0.54$, as it would be also in real observations. Due to the small number of clusters, the reliable fit subsample with $\hat{\chi}^2 < 10$ having seven clusters in total, was not applied. Therefore, in the analysis, all the 32 clusters are taken into account.

The A10 and P13 pressure profiles do not quite model the data, differently from the full sample where those two profiles are acceptable even at high redshifts (Fig. 4, Table 3), for the NIKA2 median pressure profile, see Fig. A1.

The biases from the NIKA2 subsample are compatible with the ones from the full sample, even though they are slightly higher, especially for the relaxed clusters (see Appendix A). The same problem reflects on the corrected bias, their median values are shown in the second row of Fig. A3.

These differences may be due to the selection strategy of the NIKA2 twin sample, in particular to the number of high-mass clusters, larger in proportion than the MUSIC one (for an extended analysis, see Appendix A). Anyway, based on these analyses, we can say that the NIKA2 twin sample has compatible biases within $1\sigma$ with the MUSIC sample, even if it has a larger percentage number of high-mass clusters.

# 8 CONCLUSION

The study of the bias in recovering the mass of galaxy clusters under the HE assumption is crucial for the understanding of structures formation and cosmology. In this work, the 3D ICM radial profiles of a synthetic set of almost 260 clusters from the MUSIC hydrodynamical simulations were used. The analysis has been applied on three different simulation flavours: NR, CSF, and AGN. The main difference among them is the absence of radiative processes in NR, while they are included in both CSF and AGN, where also stellar and AGN feedbacks are considered, respectively. We considered seven different redshifts in the range $0.0 \le z \le 0.82$. Moreover a MUSIC subsample, the NIKA2 twin sample, was analysed to have some hints for the NIKA2 LPSZ.

A preliminary analysis provided the clusters classification, depending on their dynamical state. In the MUSIC data set and in the NIKA2 twin sample, we have almost 50 per cent of relaxed clusters. First of all, the median pressure profile of the full sample at low redshifts ($z < 0.5$) was compared with two models describing observed local clusters from Planck (P13) and XMM-Newton (A10) in order to check whether there is consistency between simulations and the analytical models used in literature.

The SZ and X-ray estimations of the HE mass were computed from the ICM thermodynamic radial profiles using two strategies to compute the gradients: the numerical and analytical fitting. The correlation between the HE masses and the true mass $M_{500}$ was studied. Therefore, the HE mass biases were estimated. Possible dependences with redshift and dynamical state, as defined by the relaxation parameter $\chi_{DS}$, were analysed, together with the bias radial profile. We also studied the impact of non-thermal pressure support arising from the bulk motion and turbulence of the gas. All these analyses were done for both the MUSIC full data set and the NIKA2 subsample.

The main results of this work can be summarized as follows:

(i) In MUSIC, the AGN flavour better approximates the real clusters properties. In fact the pressure profile from the AGN simulation better matches the model from the observed clusters in P13 and A10.

(ii) The SZ and X-ray HE masses have a linear dependence with the true mass $M_{500}$. The slope of the fit, corresponding to the value $(1 - b)$, gives a bias $b$ of the order of 0.2, independent of the considered redshift.

(iii) Also the HE biases at $R_{500}$ are of the order of 0.2, independently from X-ray or SZ formulations and derivative estimations and they do not depend on the redshift. These results are in agreement with other simulations. The full sample and the relaxed-cluster-only biases are always in agreement with each other; the disturbed
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The mean radial profile of $\alpha = P_{\text{nth}}/P_{\text{tot}}$ for MUSIC clusters in the AGN simulation at $z=0$. The profiles for all the clusters, the only relaxed, and the only disturbed are shown as yellow, magenta, and blue solid lines, respectively. The shaded regions represent the $1\sigma$ interval.</figcaption></figure>

dependence on the redshift range analysed in this work. Moreover, the correction on $b_X$, which already gave an estimation of the mass nearer to 0 before the correction, is more effective than on the $b_{SZ}$. This leads to zero X-ray bias for the relaxed and all cluster cases, and to an overcorrection for the disturbed class. Therefore, the dynamical state has an impact on the strength of the correction, but not on its radial profile, as stated before and shown in Fig. 10.

The median radial profiles of the two corrected biases $b_{SZ,corr}$ and $b_{X,corr}$ are presented in Fig. 11 for redshift 0 and all flavours. The increasing radial dependence of the X-ray and SZ bias was cancelled when including the non-thermal correction, as expected. However, there is still a large variation in the bias profiles around $R_{500}$ for the disturbed clusters, which was already present before the correction. This could be explained by a more intriguing effect, which does not depend on the non-thermal pressure, or the simple correction formula does not account for the extreme non-thermal pressure from the disturbed clusters. We also note that deviations from spherical symmetry as well as the presence of substructures impacting the multiphase structure of the gas can play a significant role in the outskirts. As we pointed out before, the scatter (represented by the MADs in this plot) is larger than the non-corrected bias. Moreover, the SZ bias has larger scatter than the X-ray bias, especially for the CSF and NR flavours. The NR biases have less regular profiles with respect to the other flavours, but there is not a substantial difference between them. All the different flavour profiles are in agreement.

# 7 RESULTS FROM THE NIKA2 TWIN SAMPLE

The NIKA2 twin sample (Ruppin et al. 2019b) is composed of 32 MUSIC clusters, see Section 4.2. To build the sample, we chose clusters at $z = 0.82$ that are not the progenitors of the clusters at $z = 0.54$, as it would be also in real observations. Due to the small number of clusters, the reliable fit subsample with $\hat{\chi}^2 < 10$ having seven clusters in total, was not applied. Therefore, in the analysis, all the 32 clusters are taken into account.

The A10 and P13 pressure profiles do not quite model the data, differently from the full sample where those two profiles are acceptable even at high redshifts (Fig. 4, Table 3), for the NIKA2 median pressure profile, see Fig. A1.

The biases from the NIKA2 subsample are compatible with the ones from the full sample, even though they are slightly higher, especially for the relaxed clusters (see Appendix A). The same problem reflects on the corrected bias, their median values are shown in the second row of Fig. A3.

These differences may be due to the selection strategy of the NIKA2 twin sample, in particular to the number of high-mass clusters, larger in proportion than the MUSIC one (for an extended analysis, see Appendix A). Anyway, based on these analyses, we can say that the NIKA2 twin sample has compatible biases within $1\sigma$ with the MUSIC sample, even if it has a larger percentage number of high-mass clusters.

# 8 CONCLUSION

The study of the bias in recovering the mass of galaxy clusters under the HE assumption is crucial for the understanding of structures formation and cosmology. In this work, the 3D ICM radial profiles of a synthetic set of almost 260 clusters from the MUSIC hydrodynamical simulations were used. The analysis has been applied on three different simulation flavours: NR, CSF, and AGN. The main difference among them is the absence of radiative processes in NR, while they are included in both CSF and AGN, where also stellar and AGN feedbacks are considered, respectively. We considered seven different redshifts in the range $0.0 \le z \le 0.82$. Moreover a MUSIC subsample, the NIKA2 twin sample, was analysed to have some hints for the NIKA2 LPSZ.

A preliminary analysis provided the clusters classification, depending on their dynamical state. In the MUSIC data set and in the NIKA2 twin sample, we have almost 50 per cent of relaxed clusters. First of all, the median pressure profile of the full sample at low redshifts ($z < 0.5$) was compared with two models describing observed local clusters from Planck (P13) and XMM-Newton (A10) in order to check whether there is consistency between simulations and the analytical models used in literature.

The SZ and X-ray estimations of the HE mass were computed from the ICM thermodynamic radial profiles using two strategies to compute the gradients: the numerical and analytical fitting. The correlation between the HE masses and the true mass $M_{500}$ was studied. Therefore, the HE mass biases were estimated. Possible dependences with redshift and dynamical state, as defined by the relaxation parameter $\chi_{DS}$, were analysed, together with the bias radial profile. We also studied the impact of non-thermal pressure support arising from the bulk motion and turbulence of the gas. All these analyses were done for both the MUSIC full data set and the NIKA2 subsample.

The main results of this work can be summarized as follows:

(i) In MUSIC, the AGN flavour better approximates the real clusters properties. In fact the pressure profile from the AGN simulation better matches the model from the observed clusters in P13 and A10.

(ii) The SZ and X-ray HE masses have a linear dependence with the true mass $M_{500}$. The slope of the fit, corresponding to the value $(1 - b)$, gives a bias $b$ of the order of 0.2, independent of the considered redshift.

(iii) Also the HE biases at $R_{500}$ are of the order of 0.2, independently from X-ray or SZ formulations and derivative estimations and they do not depend on the redshift. These results are in agreement with other simulations. The full sample and the relaxed-cluster-only biases are always in agreement with each other; the disturbed
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Same as Fig. 8, but for the corrected biases.</figcaption></figure>

population shows very high biases, with large errors. This can be simply explained because the HE is not fully satisfied.

(iv) The different approaches to fit the ICM radial profiles have an impact on the scatter of the biases. The numerical derivative method gives on average a larger scatter with respect to the analytical fitting one, with the condition of excluding from the analysis the clusters with a non-reliable fit.

(v) The dynamical state has an impact on the scatter of the biases. Disturbed clusters show large dispersion. We find that biases depend on the continuous relaxation parameter $\chi_{DS}$, which we use to describe the dynamical state of the clusters.

(vi) While the median radial profile of the bias for the disturbed clusters shows large fluctuations with a large scatter, the relaxed and the full sample radial biases are flatter. Interestingly, in the case of these two populations the bias profile increases radially stressing the impact of non-thermal pressure contribution in the most external regions.

(vii) It is evident that non-thermal pressure contribution has an impact especially in the cluster outskirts. Correctly adding this contribution to the hydrostatic one, the estimated mass of the cluster really gets closer to the true one. Unfortunately, the scatter is larger, mainly due to the scatter in the fraction $P_{nth}/P_{tot}$, used to estimate the correction.

(viii) The mass biases from the NIKA2 subsample are compatible with the ones from the full sample, even though they are slightly higher, especially for the relaxed clusters.

The NIKA2 twin sample will help improve the analysis methodology of the ICM properties of the NIKA2 tSZ large program. Thanks to the results of this work, the sample will also help to probe different fundamental hypothesis, like spherical symmetry and HE.

Nowadays, the main systematic uncertainty associated with cluster cosmological constraints is the uncertainty on the HE mass bias. Using the non-thermal correction to obtain mass values closer to the real ones, but with a larger scatter, is not an efficient strategy, exactly as how to choose the ICM properties radial fitting approach. On the other hand, several other effects still need to be carefully considered, as temperature dishomogeneities, deviation from spherical symmetry, and the presence of substructures, in order to give a more accurate measure of the HE mass bias.

## ACKNOWLEDGEMENTS

We thank the Referee and Stefano Andreon for their useful comments. GG, MDP, and FDL acknowledge support from Sapienza Università di Roma thanks to Progetti di Ricerca Medi 2019, prot. RM11916B7540DD8D. GY acknowledges financial support by MICIU/FEDER under project grant PGC2018-094975-C21. WC acknowledges support from the European Research Council under grant number 670193 (the COSFORM project) and from the China Manned Space Program through its Space Application System. The authors wish to thank The Red Española de Supercomputación for granting computing time in the MareNostrum Supercomputer at the BSC-CSN where the MUSIC simulations have been run. FM, LP, FR, FK, and JFMP acknowledge financial support from the French National Research Agency in the framework of the ‘Investissements d’avenir’ program (ANR-15-IDEX-02). FR acknowledges financial supports provided by NASA through SAO Award Number SV2-82023 issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. ER acknowledge financial contribution from ASI-INAF n. 2017-14-H.0. VB acknowledges support by the DFG project nr. 415510302.

## DATA AVAILABILITY

The profiles analysed in this work were produced with the MUSIC simulations (Sembolini et al. 2013). These data can be accessed following the instructions on the website https://music.ft.uam.es/. The data specifically shown in this paper will be shared upon request to the authors.

## REFERENCES

* Adam R. et al., 2018, A&A, 609, A115

* Ameglio S., Borgani S., Pierpaoli E., Dolag K., Ettori S., Morandi A., 2009, MNRAS, 394, 479

* Andreon S., 2014, A&A, 570, L10

* Andreon S., Trinchieri G., Moretti A., Wang J., 2017, A&A, 606, A25

* Ansarifard S. et al., 2020, A&A, 634, A113

* Arnaud M., Pratt G. W., Piffaretti R., Böhringer H., Croston J. H., Pointecouteau E., 2010, A&A, 517, 1 (A10)
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Same as Fig. 8, but for the corrected biases.</figcaption></figure>

population shows very high biases, with large errors. This can be simply explained because the HE is not fully satisfied.

(iv) The different approaches to fit the ICM radial profiles have an impact on the scatter of the biases. The numerical derivative method gives on average a larger scatter with respect to the analytical fitting one, with the condition of excluding from the analysis the clusters with a non-reliable fit.

(v) The dynamical state has an impact on the scatter of the biases. Disturbed clusters show large dispersion. We find that biases depend on the continuous relaxation parameter $\chi_{DS}$, which we use to describe the dynamical state of the clusters.

(vi) While the median radial profile of the bias for the disturbed clusters shows large fluctuations with a large scatter, the relaxed and the full sample radial biases are flatter. Interestingly, in the case of these two populations the bias profile increases radially stressing the impact of non-thermal pressure contribution in the most external regions.

(vii) It is evident that non-thermal pressure contribution has an impact especially in the cluster outskirts. Correctly adding this contribution to the hydrostatic one, the estimated mass of the cluster really gets closer to the true one. Unfortunately, the scatter is larger, mainly due to the scatter in the fraction $P_{nth}/P_{tot}$, used to estimate the correction.

(viii) The mass biases from the NIKA2 subsample are compatible with the ones from the full sample, even though they are slightly higher, especially for the relaxed clusters.

The NIKA2 twin sample will help improve the analysis methodology of the ICM properties of the NIKA2 tSZ large program. Thanks to the results of this work, the sample will also help to probe different fundamental hypothesis, like spherical symmetry and HE.

Nowadays, the main systematic uncertainty associated with cluster cosmological constraints is the uncertainty on the HE mass bias. Using the non-thermal correction to obtain mass values closer to the real ones, but with a larger scatter, is not an efficient strategy, exactly as how to choose the ICM properties radial fitting approach. On the other hand, several other effects still need to be carefully considered, as temperature dishomogeneities, deviation from spherical symmetry, and the presence of substructures, in order to give a more accurate measure of the HE mass bias.

## ACKNOWLEDGEMENTS

We thank the Referee and Stefano Andreon for their useful comments. GG, MDP, and FDL acknowledge support from Sapienza Università di Roma thanks to Progetti di Ricerca Medi 2019, prot. RM11916B7540DD8D. GY acknowledges financial support by MICIU/FEDER under project grant PGC2018-094975-C21. WC acknowledges support from the European Research Council under grant number 670193 (the COSFORM project) and from the China Manned Space Program through its Space Application System. The authors wish to thank The Red Española de Supercomputación for granting computing time in the MareNostrum Supercomputer at the BSC-CSN where the MUSIC simulations have been run. FM, LP, FR, FK, and JFMP acknowledge financial support from the French National Research Agency in the framework of the ‘Investissements d’avenir’ program (ANR-15-IDEX-02). FR acknowledges financial supports provided by NASA through SAO Award Number SV2-82023 issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. ER acknowledge financial contribution from ASI-INAF n. 2017-14-H.0. VB acknowledges support by the DFG project nr. 415510302.

## DATA AVAILABILITY

The profiles analysed in this work were produced with the MUSIC simulations (Sembolini et al. 2013). These data can be accessed following the instructions on the website https://music.ft.uam.es/. The data specifically shown in this paper will be shared upon request to the authors.

## REFERENCES

* Adam R. et al., 2018, A&A, 609, A115

* Ameglio S., Borgani S., Pierpaoli E., Dolag K., Ettori S., Morandi A., 2009, MNRAS, 394, 479

* Andreon S., 2014, A&A, 570, L10

* Andreon S., Trinchieri G., Moretti A., Wang J., 2017, A&A, 606, A25

* Ansarifard S. et al., 2020, A&A, 634, A113

* Arnaud M., Pratt G. W., Piffaretti R., Böhringer H., Croston J. H., Pointecouteau E., 2010, A&A, 517, 1 (A10)
Barnes D. J., Vogelsberger M., Pearce F. A., Pop A.-R., Kannan R., Cao K., Kay S. T., Hernquist L., 2020, MNRAS, preprint (arXiv:2001.11508)

Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2013, ApJ, 777, 123

Biffi V., Sembolini F., De Petris M., Valdarnini R., Yepes G., Gottlöber S., 2014, MNRAS, 439, 588

Biffi V. et al., 2016, ApJ, 827, 112

Bleem L. E. et al., 2020, ApJS, 247, 25

Boehringer H., Werner N., 2009, preprint (arXiv:0907.4277)

Borgani S. et al., 2004, MNRAS, 348, 1078

Calvo M. et al., 2016, J. Low Temp. Phys., 184, 816

Capalbo V., De Petris M., De Luca F., Cui W., Yepes G., Knebe A., Rasia E., 2020, MNRAS, preprint (arXiv:2009.04565)

Carlstrom J. E., Holder G., Reese E. D., 2002, ARA&A, 40, 643

Cavaliere A., Fusco-Femiano R., 1978, A&A, 70, 677

Cialone G., De Petris M., Sembolini F., Yepes G., Baldi A. S., Rasia E., 2018, MNRAS, 477, 139

Cui W., Power C., Borgani S., Knebe A., Lewis G. F., Murante G., Poole G. B., 2017, MNRAS, 464, 2502

Cui W. et al., 2018, MNRAS, 480, 2898

De Luca F., De Petris M., Yepes G., Cui W., Knebe A., Rasia E., 2020, MNRAS, preprint (arXiv:2011.09002)

De Petris M. et al., 2020, EPJ Web Conf., 228, 00008

D'Onghia E., Navarro J. F., 2007, MNRAS, 380, L58

Eckert D., Ettori S., Pointecouteau E., Molendi S., Paltani S., Tchernin C., 2017, Astron. Nachr., 338, 293

Ettori S., Donnarumma A., Pointecouteau E., Reiprich T. H., Giodini S., Lovisari L., Schmidt R. W., 2013, Space Sci. Rev., 177, 119

Gottlöber S., Yepes G., 2007, ApJ, 664, 117

Gupta N., Saro A., Mohr J., Dolag K., Liu J., 2017, MNRAS, 469, 3069

Haggar R., Gray M. E., Pearce F. R., Knebe A., Cui W., Mostoghiu R., Yepes G., 2020, MNRAS, 492, 6074

Hasselfield M. et al., 2013, J. Cosmol. Astropart. Phys., 2013, 008

Henson M. A., Barnes D. J., Kay S. T., McCarthy I. G., Schaye J., 2017, MNRAS, 465, 3361

Hilton M. et al., 2018, ApJS, 235, 20

Hilton M. et al., 2020, ApJS, preprint (arXiv:2009.11043)

Hudson D. S., Mittal R., Reiprich T. H., Nulsen P. E. J., Andernach H., Sarazin C. L., 2010, A&A, 513, A37

Kay S. T., Peel M. W., Short C. J., Thomas P. A., Young O. E., Battye R. A., Liddle A. R., Pearce F. R., 2012, MNRAS, 422, 1999

Kitayama T., 2014, Prog. Theor. Exp. Phys., 2014, 06B111

Klypin A., Kravtsov A. V., Bullock J. S., Primack J. R., 2001, ApJ, 554, 903

Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Komatsu E. et al., 2011, ApJS, 192, 18

Koukoufilippas N., Alonso D., Bilicki M., Peacock J. A., 2020, MNRAS, 491, 5464

Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353

Lau E. T., Kravtsov A. V., Nagai D., 2009, ApJ, 705, 1129

Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2017, MNRAS, 466, 4442

Macciò A. V., Dutton A. A., Van Den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55

Makiya R., Hikage C., Komatsu E., 2020, PASJ, 72, 26

Mantz A., Allen S. W., Rapetti D., Ebeling H., 2010, MNRAS, 406, 1759

Marriage T. A. et al., 2010, ApJ, 737, 61

Martizzi D., Agrusa H. F., 2016, MNRAS, preprint (arXiv:1608.04388)

Mayet F. et al., 2020, EPJ Web Conf., 228, 00017

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

McCarthy I. G., Schaye J., Bird S., Brun A. M. C. L., 2016, MNRAS, 465, 2936

Meneghetti M., Rasia E., Merten J., Bellagamba F., Ettori S., Mazzotta P., Dolag K., Marri S., 2010, A&A, 514, A93

Mroczkowski T. et al., 2019, Space Sci. Rev., 215, 17

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

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

Nelson K., Lau E. T., Nagai D., Rudd D. H., Yu L., 2014a, ApJ, 782, 107

Neto A. F. et al., 2007, MNRAS, 381, 1450

Pacaud F. et al., 2007, MNRAS, 382, 1289

Pearce F. A., Kay S. T., Barnes D. J., Bower R. G., Schaller M., 2020, MNRAS, 491, 1622

Perotto L. et al., 2020, A&A, 637, A71

Piffaretti R., Valdarnini R., 2008, A&A, 491, 71

Planck Collaboration V, 2014a, A&A, 550, A131 (P13)

Planck Collaboration VIII, 2011a, A&A, 536, A8

Planck Collaboration XI, 2011b, A&A, 536, A11

Planck Collaboration XX, 2014b, A&A, 571, A20

Planck Collaboration XXIV, 2016a, A&A, 594, A24

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

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

Planelles S., Borgani S., Fabjan D., Killedar M., Murante G., Granato G. L., Ragone-Figueroa C., Dolag K., 2013, MNRAS, 438, 195

Prada F., Klypin A., Cuesta A., Betancort-Rijo J., Primack J., 2012, MNRAS, 423, 3018

Pratt G., Arnaud M., Biviano A., Eckert D., Ettori S., Nagai D., Okabe N., Reiprich T. H., 2019, Space Sci. Rev., 215, 25

Pratt G., Croston J., Arnaud M., Böhringer H., 2009, A&A, 498, 361

Rasia E. et al., 2006, MNRAS, 369, 2013

Rasia E. et al., 2012, New J. Phys., 14, 055018

Rasia E. et al., 2014, ApJ, 791, 96

Reichardt C. L. et al., 2013, ApJ, 763, 127

Ruppin F., Mayet F., Macías-Pérez J. F., Perotto L., 2019a, MNRAS, 490, 784

Ruppin F. et al., 2019b, A&A, 631, A21

Salvati L., Douspis M., Aghanim N., 2018, A&A, 614, A13

Sembolini F., Yepes G., De Petris M., Gottlöber S., Lamagna L., Comis B., 2013, MNRAS, 429, 323

Shi X., Komatsu E., Nagai D., Lau E. T., 2016, MNRAS, 455, 2936

Springel V., 2005, MNRAS, 364, 1105

Springel V., 2010, MNRAS, 401, 791

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

Vazza F., Angelinelli M., Jones T. W., Eckert D., Brüggen M., Brunetti G., Gheller C., 2018, MNRAS, 481, L120

Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Speybroeck L. V., 2006, ApJ, 640, 691

Vikhlinin A. et al., 2009, ApJ, 692, 1060

Williamson R. et al., 2011, ApJ, 738, 139

## APPENDIX A: RESULTS FOR THE NIKA2 TWIN SAMPLE

The median pressure profile of the NIKA2 twin sample is plotted in Fig. A1 and the parameters are listed in Table A1, first row. In this case, the profiles from A10 and P13 do not approximate well the data, even though the subsample is taken from MUSIC. To check whether this behaviour is expected, 10 subsamples were randomly extracted from MUSIC, keeping the same features in mass, redshift, and dynamical state than the NIKA2 twin sample, see Section 4. Due to the limited number of such kind of objects in MUSIC, several clusters are in common with all the 10 subsamples. In all the cases, the gNFW profile parameters are consistent ($1\sigma$), see Table A1, second row, but different to the MUSIC ones in Table 3. Repeating the same procedure, without constraining the mass range, leads to random samples that have gNFW profile parameters compatible with the MUSIC one, an example of parameters is listed in the third row of Table A1. The difference may be due to the number of high-mass clusters in the NIKA2 twin sample, larger in proportion than the MUSIC one.

In the first row of Fig. A2, the HE masses, estimated for the NIKA2 twin sample, are shown as a function of the true cluster mass. These HE masses were estimated using the analytical fitting, for both redshifts and all simulation flavours. In this figure, we just
Barnes D. J., Vogelsberger M., Pearce F. A., Pop A.-R., Kannan R., Cao K., Kay S. T., Hernquist L., 2020, MNRAS, preprint (arXiv:2001.11508)

Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2013, ApJ, 777, 123

Biffi V., Sembolini F., De Petris M., Valdarnini R., Yepes G., Gottlöber S., 2014, MNRAS, 439, 588

Biffi V. et al., 2016, ApJ, 827, 112

Bleem L. E. et al., 2020, ApJS, 247, 25

Boehringer H., Werner N., 2009, preprint (arXiv:0907.4277)

Borgani S. et al., 2004, MNRAS, 348, 1078

Calvo M. et al., 2016, J. Low Temp. Phys., 184, 816

Capalbo V., De Petris M., De Luca F., Cui W., Yepes G., Knebe A., Rasia E., 2020, MNRAS, preprint (arXiv:2009.04565)

Carlstrom J. E., Holder G., Reese E. D., 2002, ARA&A, 40, 643

Cavaliere A., Fusco-Femiano R., 1978, A&A, 70, 677

Cialone G., De Petris M., Sembolini F., Yepes G., Baldi A. S., Rasia E., 2018, MNRAS, 477, 139

Cui W., Power C., Borgani S., Knebe A., Lewis G. F., Murante G., Poole G. B., 2017, MNRAS, 464, 2502

Cui W. et al., 2018, MNRAS, 480, 2898

De Luca F., De Petris M., Yepes G., Cui W., Knebe A., Rasia E., 2020, MNRAS, preprint (arXiv:2011.09002)

De Petris M. et al., 2020, EPJ Web Conf., 228, 00008

D'Onghia E., Navarro J. F., 2007, MNRAS, 380, L58

Eckert D., Ettori S., Pointecouteau E., Molendi S., Paltani S., Tchernin C., 2017, Astron. Nachr., 338, 293

Ettori S., Donnarumma A., Pointecouteau E., Reiprich T. H., Giodini S., Lovisari L., Schmidt R. W., 2013, Space Sci. Rev., 177, 119

Gottlöber S., Yepes G., 2007, ApJ, 664, 117

Gupta N., Saro A., Mohr J., Dolag K., Liu J., 2017, MNRAS, 469, 3069

Haggar R., Gray M. E., Pearce F. R., Knebe A., Cui W., Mostoghiu R., Yepes G., 2020, MNRAS, 492, 6074

Hasselfield M. et al., 2013, J. Cosmol. Astropart. Phys., 2013, 008

Henson M. A., Barnes D. J., Kay S. T., McCarthy I. G., Schaye J., 2017, MNRAS, 465, 3361

Hilton M. et al., 2018, ApJS, 235, 20

Hilton M. et al., 2020, ApJS, preprint (arXiv:2009.11043)

Hudson D. S., Mittal R., Reiprich T. H., Nulsen P. E. J., Andernach H., Sarazin C. L., 2010, A&A, 513, A37

Kay S. T., Peel M. W., Short C. J., Thomas P. A., Young O. E., Battye R. A., Liddle A. R., Pearce F. R., 2012, MNRAS, 422, 1999

Kitayama T., 2014, Prog. Theor. Exp. Phys., 2014, 06B111

Klypin A., Kravtsov A. V., Bullock J. S., Primack J. R., 2001, ApJ, 554, 903

Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Komatsu E. et al., 2011, ApJS, 192, 18

Koukoufilippas N., Alonso D., Bilicki M., Peacock J. A., 2020, MNRAS, 491, 5464

Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353

Lau E. T., Kravtsov A. V., Nagai D., 2009, ApJ, 705, 1129

Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2017, MNRAS, 466, 4442

Macciò A. V., Dutton A. A., Van Den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55

Makiya R., Hikage C., Komatsu E., 2020, PASJ, 72, 26

Mantz A., Allen S. W., Rapetti D., Ebeling H., 2010, MNRAS, 406, 1759

Marriage T. A. et al., 2010, ApJ, 737, 61

Martizzi D., Agrusa H. F., 2016, MNRAS, preprint (arXiv:1608.04388)

Mayet F. et al., 2020, EPJ Web Conf., 228, 00017

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

McCarthy I. G., Schaye J., Bird S., Brun A. M. C. L., 2016, MNRAS, 465, 2936

Meneghetti M., Rasia E., Merten J., Bellagamba F., Ettori S., Mazzotta P., Dolag K., Marri S., 2010, A&A, 514, A93

Mroczkowski T. et al., 2019, Space Sci. Rev., 215, 17

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

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

Nelson K., Lau E. T., Nagai D., Rudd D. H., Yu L., 2014a, ApJ, 782, 107

Neto A. F. et al., 2007, MNRAS, 381, 1450

Pacaud F. et al., 2007, MNRAS, 382, 1289

Pearce F. A., Kay S. T., Barnes D. J., Bower R. G., Schaller M., 2020, MNRAS, 491, 1622

Perotto L. et al., 2020, A&A, 637, A71

Piffaretti R., Valdarnini R., 2008, A&A, 491, 71

Planck Collaboration V, 2014a, A&A, 550, A131 (P13)

Planck Collaboration VIII, 2011a, A&A, 536, A8

Planck Collaboration XI, 2011b, A&A, 536, A11

Planck Collaboration XX, 2014b, A&A, 571, A20

Planck Collaboration XXIV, 2016a, A&A, 594, A24

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

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

Planelles S., Borgani S., Fabjan D., Killedar M., Murante G., Granato G. L., Ragone-Figueroa C., Dolag K., 2013, MNRAS, 438, 195

Prada F., Klypin A., Cuesta A., Betancort-Rijo J., Primack J., 2012, MNRAS, 423, 3018

Pratt G., Arnaud M., Biviano A., Eckert D., Ettori S., Nagai D., Okabe N., Reiprich T. H., 2019, Space Sci. Rev., 215, 25

Pratt G., Croston J., Arnaud M., Böhringer H., 2009, A&A, 498, 361

Rasia E. et al., 2006, MNRAS, 369, 2013

Rasia E. et al., 2012, New J. Phys., 14, 055018

Rasia E. et al., 2014, ApJ, 791, 96

Reichardt C. L. et al., 2013, ApJ, 763, 127

Ruppin F., Mayet F., Macías-Pérez J. F., Perotto L., 2019a, MNRAS, 490, 784

Ruppin F. et al., 2019b, A&A, 631, A21

Salvati L., Douspis M., Aghanim N., 2018, A&A, 614, A13

Sembolini F., Yepes G., De Petris M., Gottlöber S., Lamagna L., Comis B., 2013, MNRAS, 429, 323

Shi X., Komatsu E., Nagai D., Lau E. T., 2016, MNRAS, 455, 2936

Springel V., 2005, MNRAS, 364, 1105

Springel V., 2010, MNRAS, 401, 791

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

Vazza F., Angelinelli M., Jones T. W., Eckert D., Brüggen M., Brunetti G., Gheller C., 2018, MNRAS, 481, L120

Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Speybroeck L. V., 2006, ApJ, 640, 691

Vikhlinin A. et al., 2009, ApJ, 692, 1060

Williamson R. et al., 2011, ApJ, 738, 139

## APPENDIX A: RESULTS FOR THE NIKA2 TWIN SAMPLE

The median pressure profile of the NIKA2 twin sample is plotted in Fig. A1 and the parameters are listed in Table A1, first row. In this case, the profiles from A10 and P13 do not approximate well the data, even though the subsample is taken from MUSIC. To check whether this behaviour is expected, 10 subsamples were randomly extracted from MUSIC, keeping the same features in mass, redshift, and dynamical state than the NIKA2 twin sample, see Section 4. Due to the limited number of such kind of objects in MUSIC, several clusters are in common with all the 10 subsamples. In all the cases, the gNFW profile parameters are consistent ($1\sigma$), see Table A1, second row, but different to the MUSIC ones in Table 3. Repeating the same procedure, without constraining the mass range, leads to random samples that have gNFW profile parameters compatible with the MUSIC one, an example of parameters is listed in the third row of Table A1. The difference may be due to the number of high-mass clusters in the NIKA2 twin sample, larger in proportion than the MUSIC one.

In the first row of Fig. A2, the HE masses, estimated for the NIKA2 twin sample, are shown as a function of the true cluster mass. These HE masses were estimated using the analytical fitting, for both redshifts and all simulation flavours. In this figure, we just
<figure><img src="image_13.png" /><figcaption>**Figure A1.** As Fig. 3, the median fit to the universal pressure profile for the NIKA2 twin sample clusters (black dots). The bottom panel shows the deviation between the median profiles for NIKA2 sample and the P13 and A10 profiles.</figcaption></figure>

show the results for the AGN simulations, both at $z = 0.54$, (open symbols), and $z = 0.82$, (solid symbols). Also in this case a fit of the type $M_{HE} = aM_{500}$ was performed, the slopes $a$ are listed in Table A2. The slopes and their errors on the NIKA2 twin sample fits are higher with respect to the full MUSIC sample, (see Table 4), due to the difference in the mass distributions. In the second row of Fig. A2, the ratio $M_{HE}/M_{500}$ is represented as a function of $M_{500}$, also in this case, as in the MUSIC sample, we do not find any bias dependence on the mass.

The HE mass biases of the NIKA2 twin sample are represented in the first row of Fig. A3 along the redshifts, using the analytical fitting method. They are slightly higher than the correspondent MUSIC cases in Fig. 9, at redshifts 0.54 and 0.82, but still in agreement with them. Contrary to MUSIC full sample, here the disturbed clusters show a bias close to 0 at $z = 0.54$. Large errors make these biases

<figure><img src="image_14.png" /><figcaption>**Figure A2.** Same as Fig. 6, but for the NIKA2 subsample. In the first row, the clusters at $z = 0.54$ are represented as open symbols, while those at $z = 0.82$ as solid symbols.</figcaption></figure>

still compatible with the whole sample. We have to report that this unexpected behaviour could be also due to having not excluded the clusters with non-reliable fits. However, the disturbed clusters show the same behaviour using the numerical derivative method.

The same problem reflects on the corrected bias. The correction on the biases has been applied to this sample as well, their median values are shown in the second row of Fig. A3. We see that, at redshift 0.54, the non-thermal correction works well, giving bias values close to 0, compatible with the MUSIC sample. Instead, at redshift 0.82 the corrected bias is slightly lower than 0, although still compatible with it within $2\sigma$. At $z = 0.54$, the correction does not work for the case of disturbed clusters, but it is still in agreement with the MUSIC value.

<table><caption>**Table A1.** Best-fitting parameters and errors to the gNFW pressure profile for the NIKA2 twin sample (first row). The second row shows the same information but for a random selection of MUSIC clusters with the same mass range, redshift, and dynamical state than the NIKA2 sample, while in the third row, the random selection is done without any mass constrain.</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>NIKA2 twin sample</td><td>10.21 ± 2.55</td><td>0.19 ± 0.10</td><td>0.78 ± 0.06</td><td>12.13 ± 3.10</td><td>0.20 ± 0.06</td></tr><tr><td>Random sample</td><td>16.89 ± 5.37</td><td>0.18 ± 0.09</td><td>0.77 ± 0.06</td><td>12.92 ± 3.10</td><td>0.09 ± 0.06</td></tr><tr><td>Random sample (no mass cut)</td><td>3.79 ± 0.35</td><td>2.04 ± 0.07</td><td>1.82 ± 0.12</td><td>3.71 ± 0.07</td><td>0.56 ± 0.03</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure A1.** As Fig. 3, the median fit to the universal pressure profile for the NIKA2 twin sample clusters (black dots). The bottom panel shows the deviation between the median profiles for NIKA2 sample and the P13 and A10 profiles.</figcaption></figure>

show the results for the AGN simulations, both at $z = 0.54$, (open symbols), and $z = 0.82$, (solid symbols). Also in this case a fit of the type $M_{HE} = aM_{500}$ was performed, the slopes $a$ are listed in Table A2. The slopes and their errors on the NIKA2 twin sample fits are higher with respect to the full MUSIC sample, (see Table 4), due to the difference in the mass distributions. In the second row of Fig. A2, the ratio $M_{HE}/M_{500}$ is represented as a function of $M_{500}$, also in this case, as in the MUSIC sample, we do not find any bias dependence on the mass.

The HE mass biases of the NIKA2 twin sample are represented in the first row of Fig. A3 along the redshifts, using the analytical fitting method. They are slightly higher than the correspondent MUSIC cases in Fig. 9, at redshifts 0.54 and 0.82, but still in agreement with them. Contrary to MUSIC full sample, here the disturbed clusters show a bias close to 0 at $z = 0.54$. Large errors make these biases

<figure><img src="image_15.png" /><figcaption>**Figure A2.** Same as Fig. 6, but for the NIKA2 subsample. In the first row, the clusters at $z = 0.54$ are represented as open symbols, while those at $z = 0.82$ as solid symbols.</figcaption></figure>

still compatible with the whole sample. We have to report that this unexpected behaviour could be also due to having not excluded the clusters with non-reliable fits. However, the disturbed clusters show the same behaviour using the numerical derivative method.

The same problem reflects on the corrected bias. The correction on the biases has been applied to this sample as well, their median values are shown in the second row of Fig. A3. We see that, at redshift 0.54, the non-thermal correction works well, giving bias values close to 0, compatible with the MUSIC sample. Instead, at redshift 0.82 the corrected bias is slightly lower than 0, although still compatible with it within $2\sigma$. At $z = 0.54$, the correction does not work for the case of disturbed clusters, but it is still in agreement with the MUSIC value.

<table><caption>**Table A1.** Best-fitting parameters and errors to the gNFW pressure profile for the NIKA2 twin sample (first row). The second row shows the same information but for a random selection of MUSIC clusters with the same mass range, redshift, and dynamical state than the NIKA2 sample, while in the third row, the random selection is done without any mass constrain.</caption><thead><tr><th></th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>NIKA2 twin sample</td><td>10.21 ± 2.55</td><td>0.19 ± 0.10</td><td>0.78 ± 0.06</td><td>12.13 ± 3.10</td><td>0.20 ± 0.06</td></tr><tr><td>Random sample</td><td>16.89 ± 5.37</td><td>0.18 ± 0.09</td><td>0.77 ± 0.06</td><td>12.92 ± 3.10</td><td>0.09 ± 0.06</td></tr><tr><td>Random sample (no mass cut)</td><td>3.79 ± 0.35</td><td>2.04 ± 0.07</td><td>1.82 ± 0.12</td><td>3.71 ± 0.07</td><td>0.56 ± 0.03</td></tr></tbody></table>
<table><caption>**Table A2.** The best-fitting values and their corresponding errors to the parameter $a = 1 - b$ in the $M_{HE} = aM_{500}$ relation for the NIKA2 twin sample in the two redshift bins and for all the simulation flavours.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Slope a<sub>SZ</sub></th><th colspan="3">Slope a<sub>X</sub></th></tr><tr><th>AGN</th><th>CSF</th><th>NR</th><th>AGN</th><th>CSF</th><th>NR</th></tr></thead><tbody><tr><td>0.54+0.82</td><td>0.783 ± 0.028</td><td>0.948 ± 0.071</td><td>0.775 ± 0.034</td><td>0.970 ± 0.037</td><td>0.952 ± 0.030</td><td>0.972 ± 0.050</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure A3.** The redshift evolution of the median values of the SZ and X-ray bias (first row) and the corrected ones, (second row), for the NIKA2 twin sample, using the AGN flavour and the analytical fitting method. The relaxed, disturbed, and all clusters cases are represented in magenta diamonds, blue squares, and yellow triangles, respectively. The error bars represent the 16th and 84th percentile.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure B1.** Histograms of the $b_{SZ}$ and $b_X$ values corresponding to AGN simulations of the complete MUSIC sample, using the analytical fitting method, for the whole MUSIC sample, without discriminating for the goodness of the fit $\hat{\chi}^2$. In dark cyan (left figure) and light green (right figure), there are the biases of all the clusters, in magenta, for only relaxed and in blue for only disturbed clusters. In the legend, the values of the skewness and the kurtosis are reported.</figcaption></figure>

## APPENDIX B: MASS BIAS DISTRIBUTIONS

Examples of bias distribution are shown in Fig. B1 for the analytical fitting method, without discriminating for the goodness of the fits $\hat{\chi}^2$, for the AGN simulation at $z = 0$. In Table B1, we report the means and the standard deviations from the data and from a Gaussian fit for both methods. The distributions are not Gaussian, as evident from the skewness and kurtosis values, reported in the legends of Fig. B1. The Anderson test also confirms that the distributions are not Gaussian and thus we will not comment on the mean and standard deviation. Only in a few cases the distributions are actually Gaussian. In the Table the bootstrap errors, $err_b$, are listed also. For all cases they are of the order of $10^{-2}$.

We can say that the two methods to compute the derivative lead to very similar and compatible results, nevertheless it is worth noticing that the biases from the numerical derivative method show a broader distribution. Computing the mass bias through the X-ray and SZ equations based on different ICM properties also leads to similar results on the HE bias. However, in the case of the analytical fitting method, $b_{SZ}$ and $b_X$ have both negative skewness for all the redshifts, meaning that they have a more pronounced tail on the left side of the distribution, i.e. towards the 0 bias, as found also in Ansarifard et al. (2020). They study the distributions of the $b_{SZ}$ and $b_X$, estimated using the gas density $\rho_g$ instead of $\mu m_p n_e$ in equations (1) and (2), finding that $b_{SZ}$ shows a larger scatter and that $b_X$ is more symmetric than $b_{SZ}$. These results disagree with what we found here. We conclude that the differences are due to the radial gas density profiles considered to estimate $M_{HE}$. In fact, if we use $\rho_g$, we have the same distribution features as in Ansarifard et al. (2020).
<table><caption>**Table A2.** The best-fitting values and their corresponding errors to the parameter $a = 1 - b$ in the $M_{HE} = aM_{500}$ relation for the NIKA2 twin sample in the two redshift bins and for all the simulation flavours.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Slope a<sub>SZ</sub></th><th colspan="3">Slope a<sub>X</sub></th></tr><tr><th>AGN</th><th>CSF</th><th>NR</th><th>AGN</th><th>CSF</th><th>NR</th></tr></thead><tbody><tr><td>0.54+0.82</td><td>0.783 ± 0.028</td><td>0.948 ± 0.071</td><td>0.775 ± 0.034</td><td>0.970 ± 0.037</td><td>0.952 ± 0.030</td><td>0.972 ± 0.050</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Figure A3.** The redshift evolution of the median values of the SZ and X-ray bias (first row) and the corrected ones, (second row), for the NIKA2 twin sample, using the AGN flavour and the analytical fitting method. The relaxed, disturbed, and all clusters cases are represented in magenta diamonds, blue squares, and yellow triangles, respectively. The error bars represent the 16th and 84th percentile.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure B1.** Histograms of the $b_{SZ}$ and $b_X$ values corresponding to AGN simulations of the complete MUSIC sample, using the analytical fitting method, for the whole MUSIC sample, without discriminating for the goodness of the fit $\hat{\chi}^2$. In dark cyan (left figure) and light green (right figure), there are the biases of all the clusters, in magenta, for only relaxed and in blue for only disturbed clusters. In the legend, the values of the skewness and the kurtosis are reported.</figcaption></figure>

## APPENDIX B: MASS BIAS DISTRIBUTIONS

Examples of bias distribution are shown in Fig. B1 for the analytical fitting method, without discriminating for the goodness of the fits $\hat{\chi}^2$, for the AGN simulation at $z = 0$. In Table B1, we report the means and the standard deviations from the data and from a Gaussian fit for both methods. The distributions are not Gaussian, as evident from the skewness and kurtosis values, reported in the legends of Fig. B1. The Anderson test also confirms that the distributions are not Gaussian and thus we will not comment on the mean and standard deviation. Only in a few cases the distributions are actually Gaussian. In the Table the bootstrap errors, $err_b$, are listed also. For all cases they are of the order of $10^{-2}$.

We can say that the two methods to compute the derivative lead to very similar and compatible results, nevertheless it is worth noticing that the biases from the numerical derivative method show a broader distribution. Computing the mass bias through the X-ray and SZ equations based on different ICM properties also leads to similar results on the HE bias. However, in the case of the analytical fitting method, $b_{SZ}$ and $b_X$ have both negative skewness for all the redshifts, meaning that they have a more pronounced tail on the left side of the distribution, i.e. towards the 0 bias, as found also in Ansarifard et al. (2020). They study the distributions of the $b_{SZ}$ and $b_X$, estimated using the gas density $\rho_g$ instead of $\mu m_p n_e$ in equations (1) and (2), finding that $b_{SZ}$ shows a larger scatter and that $b_X$ is more symmetric than $b_{SZ}$. These results disagree with what we found here. We conclude that the differences are due to the radial gas density profiles considered to estimate $M_{HE}$. In fact, if we use $\rho_g$, we have the same distribution features as in Ansarifard et al. (2020).