# Euclid preparation - XXXIX. The effect of baryons on the Halo Mass Function

T Castro, S Borgani, M Costanzi, J Dakin, K Dolag, A Fumagalli, A Ragagnin, A Saro, A.M.C Le Brun, Nabila Aghanim, et al.

► To cite this version:

T Castro, S Borgani, M Costanzi, J Dakin, K Dolag, et al., Euclid preparation - XXXIX. The effect of baryons on the Halo Mass Function. Astron.Astrophys., 2024, 685, pp.A109. 10.1051/0004-6361/202348388. [arXiv:2004.03617]

HAL Id: hal-04273972

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

Submitted on 18 Oct 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
Euclid preparation

XXXIX. The effect of baryons on the halo mass function

Euclid Collaboration: T. Castro1,2,3,4, S. Borgani5,3,1,2, M. Costanzi1,5, J. Dakin6, K. Dolag7, A. Fumagalli5,1,2,3, A. Ragagnin8,3,9, A. Saro5,3,1,2,4, A. M. C. Le Brun10, N. Aghanim11, A. Amara12, S. Andreon13, N. Auricchio8, M. Baldi14,8,15, S. Bardelli8, C. Bodendorf16, D. Bonino17, E. Branchini18,19, M. Brescia20,21,22, J. Brinchmann23, S. Camera24,25,17, V. Capobianco17, C. Carbone26, J. Carretero27,28, S. Casas29, M. Castellano30, S. Cavuoti21,22, A. Cimatti31, G. Congedo32, C. J. Conselice33, L. Conversi34,35, Y. Copin36, L. Corcione17, F. Courbin37, H. M. Courtois38, M. Cropper39, A. Da Silva40,41, H. Degaudenzi42, A. M. Di Giorgio43, J. Dinis41,40, F. Dubath42, C. A. J. Duncan44,33, X. Dupac35, M. Farina43, S. Farrens45, S. Ferriol36, M. Frailis1, E. Franceschi8, M. Fumana26, S. Galeotta1, B. Gillis32, C. Giocoli8,15, A. Grazian46, F. Grupp16,47, S. V. H. Haugan48, W. Holmes49, F. Hormuth50, A. Hornstrup51,52, K. Jahnke53, E. Keihänen54, S. Kermiche55, A. Kiessling49, M. Kilbinger56, B. Kubik36, M. Kunz57, H. Kurki-Suonio58,59, S. Ligori17, P. B. Lilje48, V. Lindholm58,59, I. Lloro60, E. Maiorano8, O. Mansutti1, O. Marggraf61, K. Markovic49, N. Martinet62, F. Marulli9,8,15, R. Massey63, S. Maurogordato64, E. Medinaceli8, M. Meneghetti8,15, E. Merlin30, G. Meylan37, M. Moresco9,8, L. Moscardini9,8,15, E. Munari1, S.-M. Niemi65, C. Padilla27, S. Paltani42, F. Pasian1, V. Pettorino66, S. Pires45, G. Polenta67, M. Poncet68, L. A. Popa69, L. Pozzetti8, F. Raison16, R. Rebolo70,71, A. Renzi72,73, J. Rhodes49, G. Riccio21, E. Romelli1, M. Roncarelli8, R. Saglia7,16, D. Sapone74, B. Sartoris7,1, P. Schneider61, T. Schrabback75, A. Secroun55, G. Seidel53, S. Serrano76,77,78, C. Sirignano72,73, G. Sirri15, L. Stanco73, J.-L. Starck56, P. Tallada-Crespi79,28, A. N. Taylor32, I. Tereno40,80, R. Toledo-Moreo81, F. Torradeflot28,79, I. Tutusaus82, E. A. Valentijn83, L. Valenziano8,84, T. Vassallo7,1, A. Veropalumbo13,19, Y. Wang85, J. Weller7,16, A. Zacchei1,3, G. Zamorani8, J. Zoubian55, E. Zucca8, A. Biviano1,3, E. Bozzo42, C. Cerna86,87, C. Colodro-Conde70, D. Di Ferdinando15, N. Mauri31,15, C. Neissner27,28, Z. Sakr88,82,89, V. Scottez90,91, M. Tenti15, M. Viel3,1,92,2,4, M. Wiesmann48, Y. Akrami93,94, S. Anselmi72,73, C. Baccigalupi92,1,2,3, M. Ballardini95,96,8, A. S. Borlaff97,98,99, S. Bruton100, C. Burigana101,84, R. Cabanac82, A. Cappi8,64, C. S. Carvalho80, G. Castignani9,8, G. Cañas-Herrera65,102, K. C. Chambers103, A. R. Cooray104, J. Coupon42, O. Cucciati8, A. Díaz-Sánchez105, S. Davini19, S. de la Torre62, G. De Lucia1, G. Desprez106, S. Di Domizio18,19, H. Dole11, S. Escoffier55, I. Ferrero48, F. Finelli8,84, L. Gabarra73,72, K. Ganga107, J. Garcia-Bellido93, F. Giacomini15, G. Gozalias108,58, H. Hildebrandt109, S. Ilić110,68,82, A. Jimanez Munñoz111, J. J. E. Kajava112,113, V. Kansal114,115, C. C. Kirkpatrick54, L. Legrand57, A. Loureiro116,117, J. Macias-Perez111, M. Magliocchetti43, G. Mainetti118, R. Maoli119,30, M. Martinelli30,120, C. J. A. P. Martins121,23, S. Matthew32, M. Maturi88,122, L. Maurin11, R. B. Metcalf9,8, M. Migliaccio123,124, P. Monaco5,1,2,3, G. Morgante8, S. Nadathur12, L. Patrizii15, A. Pezzotta16, V. Popa69, C. Porciani61, D. Potter6, M. Pöntinen58, P. Reimberg90, P.-F. Rocci11, A. G. Sánchez16, J. Schaye125, A. Schneider6, E. Sefusatti1,3,2, M. Sereno8,15, P. Simon61, A. Spurio Mancini39, J. Stadel6, S. A. Stanford126, J. Steinwagner16, G. Testera19, M. Tewes61, R. Teyssier127, S. Toft52,128, S. Tosi18,19,13, A. Troja72,73, M. Tucci42, J. Valiviita58,59, and D. Vergani8

(Affiliations can be found after the references)

Received 25 October 2023 / Accepted 31 January 2024

ABSTRACT

The *Euclid* photometric survey of galaxy clusters stands as a powerful cosmological tool, with the capacity to significantly propel our understanding of the Universe. Despite being subdominant to dark matter and dark energy, the baryonic component of our Universe holds substantial influence over the structure and mass of galaxy clusters. This paper presents a novel model that can be used to precisely quantify the impact of baryons on the virial halo masses of galaxy clusters using the baryon fraction within a cluster as a proxy for their effect. Constructed on the premise of quasi-adiabaticity, the model includes two parameters, which are calibrated using non-radiative cosmological hydrodynamical simulations, and a single large-scale simulation from the Magneticum set, which includes the physical processes driving galaxy formation. As a main result of our analysis, we demonstrate that this model delivers a remarkable 1% relative accuracy in determining the virial dark matter-only equivalent mass of galaxy clusters starting from the corresponding total cluster mass and baryon fraction measured in hydrodynamical simulations. Furthermore, we demonstrate that this result is robust against changes in cosmological parameters and against variation of the numerical implementation of the subresolution physical processes included in the simulations. Our work substantiates previous claims regarding the impact of baryons on cluster cosmology studies. In particular, we show how neglecting these effects would lead to biased cosmological constraints for a *Euclid*-like cluster abundance analysis. Importantly, we demonstrate that uncertainties associated with our model arising from baryonic corrections to cluster masses are subdominant when compared to the precision with which mass-observable (i.e. richness) relations will be calibrated using *Euclid* and to our current understanding of the baryon fraction within galaxy clusters.

**Key words.** galaxies: clusters: general – cosmology: theory – large-scale structure of Universe
# 1. Introduction

In the Lambda Cold Dark Matter (ACDM) scenario, structures are formed hierarchically, with larger systems emerging from the merger of smaller ones. As the largest, most massive, virialized objects in the Universe, galaxy clusters are at the pinnacle of this hierarchy and provide competitive cosmological probes of the geometry of our Universe and of the growth of density perturbations through measurements of their cosmic abundance and clustering (see e.g. Allen et al. 2011; Kravtsov & Borgani 2012; Fumagalli et al. 2023).

While structure formation is gravitationally dominated by dark matter, astrophysical effects associated with the baryonic component are known to significantly affect the clusters and cluster galaxy population (see e.g. McDonald et al. 2012; Webb et al. 2015; Ellien et al. 2019; Schellenberger et al. 2019; Yuan et al. 2020; Debackere et al. 2021). In this context, advanced cosmological simulations (e.g. Borgani & Kravtsov 2011) provide the best tools for studying the formation of galaxy clusters starting from primordial density perturbations.

The primary cosmological probe from cluster surveys is derived from number-count experiments (Borgani et al. 2001; Holder et al. 2001; Rozo et al. 2010; Hasselfield et al. 2013; Planck Collaboration XX 2014; Bocquet et al. 2015, 2019; Mantz et al. 2015; Planck Collaboration XXIV 2016; Abbott et al. 2020; Costanzi et al. 2021; Lesci et al. 2022), which rely on the strong dependence of both volumes and the density of halos as a function of mass and redshift on cosmological parameters. The number density of halos as a function of mass and redshift is called the halo mass function (HMF).

The prediction of the HMF relies deeply on the results of cosmological simulations (see e.g. Tinker et al. 2008; Watson et al. 2013; Bocquet et al. 2016, 2020; Despali et al. 2016; Castro et al. 2020; Euclid Collaboration 2023). As simulations involving full hydrodynamical calculations are more expensive than purely gravitational N-body simulations, and the cost can even be prohibitive, the conventional method is to characterise the HMF using the latter and model how baryonic physics alters the mass of halos in post-processing (Schneider & Teyssier 2015; Aricò et al. 2021). For the sake of brevity, we refer to these hydrodynamical simulations simply as ‘hydro’ and their dark-matter(DM)-only N-body counterparts as ‘dmo’. As mentioned, despite being sub-dominant, the baryonic component has a sizeable impact on the detailed properties of the large-scale structure of our Universe (Cui et al. 2014; Velliscig et al. 2014; Bocquet et al. 2016; Castro et al. 2018, 2020; Schaye et al. 2023).

Feedback processes related to supernovae (SNe) and to active galactic nuclei (AGN) originate on small scales that cannot be explicitly resolved in simulations covering large cosmological volumes. As a result, it is impossible to explicitly simulate these processes within such large volumes from first principles; instead, subresolution models must be used (e.g. Springel & Hernquist 2003; Hirschmann et al. 2014; Vogelsberger et al. 2014, 2020; Schaye et al. 2015; Crain et al. 2015; McCarthy et al. 2017). However, these processes affect the distribution of baryons and cosmic structures on scales well resolved by cosmological simulations and probed by observations (see e.g. van Daalen et al. 2011).

The strategy of describing such processes in simulations through subresolution phenomenological models has offered invaluable insights for quantifying the influence of baryonic effects on large-scale structure (see e.g. Teyssier et al. 2011;

van Daalen et al. 2011, 2020; Martizzi et al. 2013; Sawala et al. 2016; Bullock & Boylan-Kolchin 2017; Schaye et al. 2023). However, due to our incomplete understanding of the associated subresolution physical processes, several of the parameters used in their implementation require a calibration to reproduce some key observables. This calibration is usually carried out for one or a few specific observables, and is intrinsically resolution dependent (see Schaye et al. 2015, for discussion). Furthermore, due to the complex interaction and the degeneracy between highly non-linear processes, different choices for such parameters can equally reproduce the target observables.

It is established that at the scales relevant for galaxy clusters, baryonic feedback cannot disrupt structures. However, various astrophysical processes associated with the baryonic component can still alter the halo mass. Specifically, these processes include radiative cooling and star formation, both of which can induce adiabatic contraction of the halo. Moreover, sudden displacement of gas caused by impulsive AGN feedback can lead to expansion of the halo. These dynamic processes significantly modify the halo mass when compared to the same halo described within a dmo simulation (see e.g. Velliscig et al. 2014; Castro et al. 2020). In this light, a comprehensive understanding of such influences is indispensable in order to accurately model and interpret hydro simulations.

In the present paper, we model the impact of baryonic effects on the halo virial mass. Our model assumes that the processes involved are quasi-adiabatic; that is, that they happen in a steady, non-abrupt manner. Although our model is calibrated on a set of non-radiative hydro simulations and on a single realization of a full-physics simulation, we demonstrate that its performance is robust against the change of parameters controlling the subresolution physics and our choice between independent implementations of subresolution processes.

The large statistics of clusters that will be provided by the Euclid (Laureijs et al. 2011) cluster survey require a high-precision calibration of the HMF, so that uncertainties on this theory-predicted quantity do not represent a limiting factor in the precision achievable by the cluster survey (see Euclid Collaboration 2019). In Euclid Collaboration (2023), we used a large set of N-body simulations to derive an analytical, cosmology-dependent expression for the HMF that meets the required precision of ≈1%. In the present paper, we extend the results of our previous work to account for baryonic effects. Our goal is therefore to reach a 1% accuracy on the baryonic impact description in order to match the Euclid requirements and allow optimum exploitation of the cosmological constraining power of the Euclid cluster survey (Sartoris et al. 2016).

This paper is organized as follows: we outline the methodology used in this paper in Sect. 2. In Sect. 3, we present our findings from the simulations and the performance of the baryonic correction model. In Sect. 4, we assess the impact of our model in a forecast Euclid cluster-count analysis. Some final remarks are provided in Sect. 5. Lastly, in Appendix A, we study the baryonic content in different simulations.

# 2. Methodology

## 2.1. Simulations

In this section, we describe the simulations used to calibrate the baryonic effects on halo masses, and the preparation of the halo catalogues used to compare results from hydro and dmo simulations.
<table><caption>**Table 1.** Set of boxes from Magneticum simulations and halo-selection parameters used in this work.</caption><thead><tr><th rowspan="2">Box</th><th rowspan="2">L<sub>box</sub><br>(Mpc)</th><th rowspan="2">m<sub>DM</sub><br>(M<sub>☉</sub>)</th><th rowspan="2">M<sub>halo,min</sub><br>(M<sub>☉</sub>)</th><th colspan="3">Physics</th></tr><tr><th>Hydro</th><th>dmo</th><th>Non-radiative</th></tr></thead><tbody><tr><td>2</td><td>500</td><td>9.8 × 10<sup>8</sup></td><td>4.3 × 10<sup>13</sup></td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>2b</td><td>909</td><td>9.8 × 10<sup>8</sup></td><td>4.3 × 10<sup>13</sup></td><td>✓<sup>(1)</sup></td><td>✓</td><td>χ</td></tr><tr><td>1a</td><td>1272</td><td>1.9 × 10<sup>10</sup></td><td>4.3 × 10<sup>14</sup></td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

**Notes.** Column 1: name of the boxes; Col. 2: size of the simulation box; Col. 3: mass of the DM particle in each box; Col. 4: minimum mass of the halos considered in our analysis; Col. 5: the different physical models used to run the simulations. <sup>(1)</sup>Box 2b ran only until $z = 0.2$.

### 2.1.1. The Magneticum set

The Magneticum<fnref n="1" /> simulations were carried out using the TreePM+SPH code P-Gadget3, which is a more efficient variation of Gadget-2 (Springel et al. 2001b; Springel 2005). The SPH solver uses the revised implementation of Beck et al. (2016), which overcomes a number of limitations of the traditional SPH solvers. The Magneticum hydro simulations include treatment of radiative cooling, heating by a uniform, evolving UV background, star formation and stellar feedback following Springel et al. (2005), and the description of stellar evolution and chemical enrichment processes by Tornatore et al. (2007). According to the latter, 11 chemical elements are followed, once produced by AGB stars, Type-Ia SN, and Type-II SN (H; He; C; N; O; Ne; Mg; Si; S; Ca; Fe). Following Wiersma et al. (2009), metallicity-dependent cooling is implemented using cooling tables generated by the freely accessible CLOUDY photo-ionization algorithm (Ferland et al. 1998). Supermassive black holes, which are hosted at the centre of galaxies, are described as sink particles whose mass increases by gas accretion and merging with other BHs (Springel et al. 2005; Di Matteo et al. 2008). Accretion onto these black holes is governed by the Bondi-Hoyle-Lyttleton formula (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), capped by the Eddington rate, and involves both quasar- and radio-mode feedback regimes. Further details on the specific model of gas accretion and the AGN feedback are provided in Hirschmann et al. (2014). The Magneticum set includes a dmo counterpart for every hydro simulation. Additionally, a select number of these hydro simulations are paired with non-radiative simulations. The non-radiative simulations include gas hydrodynamics but no other subresolution physics.

In Table 1, we report the relevant parameters regarding box size and mass resolution of the subset of Magneticum simulations we use as reference for our analysis. While the reference cosmology of the Magneticum simulations is WMAP7 (Komatsu et al. 2011), for the largest box here considered, Box 1a, the Magneticum suite also covers different cosmologies, as presented in Table 2 (see also Singh et al. 2020). Although cosmological parameters are varied, all these simulations keep the relevant parameters describing the subresolution models fixed. In addition, we also include four simulations that either vary the AGN feedback efficiency by ±33 percent, or assume a wind velocity of 500 or 800 km s⁻<fnref n="1" />, instead of the fiducial value of 350 km s⁻<fnref n="1" />.

<table><caption>**Table 2.** Set of cosmologies covered by Magneticum Box 1a simulations used in this work.</caption><thead><tr><th rowspan="2">Name</th><th rowspan="2">Ω<sub>m</sub></th><th rowspan="2">σ<sub>8</sub></th><th rowspan="2">h</th><th rowspan="2">f<sub>b</sub></th><th colspan="3">Physics</th></tr><tr><th>Hydro</th><th>dmo</th><th>Non-radiative</th></tr></thead><tbody><tr><td>C1</td><td>0.153</td><td>0.614</td><td>0.666</td><td>0.267</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>C2</td><td>0.189</td><td>0.697</td><td>0.703</td><td>0.241</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C3</td><td>0.204</td><td>0.739</td><td>0.689</td><td>0.214</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C4</td><td>0.200</td><td>0.850</td><td>0.730</td><td>0.208</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C5</td><td>0.222</td><td>0.793</td><td>0.676</td><td>0.190</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C6</td><td>0.232</td><td>0.687</td><td>0.670</td><td>0.178</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C7</td><td>0.268</td><td>0.721</td><td>0.699</td><td>0.168</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C8<sup>(1)</sup></td><td>0.272</td><td>0.809</td><td>0.704</td><td>0.168</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>C9</td><td>0.304</td><td>0.886</td><td>0.740</td><td>0.166</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C10</td><td>0.301</td><td>0.834</td><td>0.708</td><td>0.153</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C11</td><td>0.342</td><td>0.834</td><td>0.708</td><td>0.135</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C12</td><td>0.363</td><td>0.884</td><td>0.729</td><td>0.135</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C13</td><td>0.400</td><td>0.650</td><td>0.675</td><td>0.121</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C14</td><td>0.406</td><td>0.867</td><td>0.712</td><td>0.115</td><td>✓</td><td>✓</td><td>χ</td></tr><tr><td>C15</td><td>0.428</td><td>0.830</td><td>0.732</td><td>0.115</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

**Notes.** <sup>(1)</sup> The reference WMAP7 cosmology.

### 2.1.2. Illustris TNG

The Illustris TNG<fnref n="2" /> simulations were run using the AREPO code (Springel 2010), which has a gravity solver similar to the P-Gadget3 one, with a hydro solver based on solving the Riemann problem on a moving mesh. This set of simulations represents an improvement – in terms of the sub-resolution models adopted for star-formation and feedback models – with respect to those adopted in the original Illustris simulations (Vogelsberger et al. 2013). It makes use of the same galaxy formation model with an improved kinetic AGN feedback model, a new parameterisation of galactic winds, and the addition of magnetic fields (see, Weinberger et al. 2017; Pillepich et al. 2018a; Pakmor et al. 2011, 2014; Pakmor & Springel 2013, respectively).

For the purpose of the analysis presented in this paper, we use the largest box available; that is, the TNG300 (Naiman et al. 2018; Nelson et al. 2018, 2019; Springel et al. 2018; Pillepich et al. 2018b) and its dmo counterpart. The simulated box size is 302.6 Mpc and has a DM particle mass of $5.9 \times 10^7 M_{\odot}$, the assumed cosmological model being that of Planck 2015 (Planck Collaboration XIII 2016). We select all halos with mass higher than $1.0 \times 10^{13} M_{\odot} h^{-1}$.

### 2.1.3. Halo catalogues

The halo catalogues used in this paper were obtained using the SUBFIND halo finder (Springel et al. 2001a, 2021), as presented in Dolag et al. (2009). This version of the code also includes gas and star particles. SUBFIND starts by establishing halo centres by running a parallel implementation of the 3D friends-of-friends (FOF; see e.g. Davis et al. 1985) algorithm and then allocating it to the position of the particle with the lowest potential. After that, it grows spheres around the centre and calculates the mass inside several radii, including a mean density $\Delta$ times the cosmic critical density. Throughout this paper, we only refer to the quantities with respect to the virial overdensity, $\Delta_{vir}$, computed for the corresponding cosmology (see Bryan & Norman 1998).

<footnote n="1">¹ http://www.magneticum.org</footnote>

<footnote n="2">² https://www.tng-project.org</footnote>
### 2.1.4. Matching algorithm

For the purpose of the analysis presented in this paper, we need to identify the halo pairs corresponding to the same object identified in a hydro simulation and in its dmo counterpart. The matched catalogues are created by detecting the halo pairs from each of the hydro and dmo runs that are the closest based on the spatial proximity of their halo centres. Pairs that had a difference in mass of greater than 50 percent were discarded. The same procedure was applied by Castro et al. (2020), who stated that the matched catalogue completeness and purity are greater than 95 percent for objects more massive than a few times $10^{13} M_{\odot}$. We independently validated the performance of the matching algorithm by comparing this position-based matching algorithm with a stricter matching algorithm based on matching objects that share more than 50 percent of the DM particles.

## 2.2. Baryonic correction on halo masses

The baryonic component can impact the spherical overdensity mass of galaxy clusters in three ways: due to adiabatic contraction, baryonic feedback, and the different dynamics between baryonic matter and DM. The adiabatic contraction is induced by radiative cooling of the diffuse gas. As the gas loses energy, its pressure support diminishes, leading to a fast collapse towards the bottom of the potential well, which in turn reacts by slightly contracting (e.g. Gnedin et al. 2004). This contraction is modulated by the rate of gas cooling.

Concerning the feedback effects, these can cause sudden gas heating, and its subsequent displacement by expansion. Consequently, the DM halo responds by slightly expanding, with the extent of the expansion depending on the intensity of feedback. While feedback associated with SNe has been shown in simulations not to efficiently counteract adiabatic contraction on the scales of galaxy clusters and groups, AGN feedback is able to displace a large amount of gas. Consequently, the overall structure of the cluster is affected as the DM component is passively dragged by the baryonic component.

Lastly, another baryonic effect is the difference in accretion dynamics between DM and gas, with baryonic matter experiencing slower accretion due to the presence of shocks. This differential accretion can also affect the overall structure and mass of the galaxy clusters.

If baryonic effects happen slowly with respect to the cluster dynamical time and orbits are spherically symmetric, the angular momentum is an adiabatic invariant of the mass enclosed by the collapsing shell, which means there is conservation of the quantity $M(R)R$, where $M(R)$ is the total mass inside the radius $R$ (see e.g. Steigman et al. 1978; Zeldovich et al. 1980; Blumenthal et al. 1986; Ryden & Gunn 1987). Therefore, for a cluster-sized halo produced in the hydrodynamic simulation with virial mass $M_{vir,hyd}$ and virial radius $R_{vir,hyd}$ we can write

$$M_{\Delta,dmo} R_{\Delta,dmo} = M_{vir,hyd} R_{vir,hyd}, \quad (1)$$

where $M_{\Delta,dmo}$ is the mass contained inside the sphere of radius $R_{\Delta,dmo}$ centred on the dmo counterpart not affected by the baryonic physics. We note that, due to the baryonic impact on the halo profile, the overdensity threshold $\Delta$ in general will differ from the virial value. Defining $f_{b,cosmic}$ as the cosmic baryon fraction and $f_{b,vir}$ as the baryon fraction within the halo virial radius, then

$$M_{\Delta,dmo} = \frac{1 - f_{b,vir}}{1 - f_{b,cosmic}} M_{vir,hyd}. \quad (2)$$

Equation (2) is derived from the consideration that the DM contribution to the spherical overdensity mass in hydro-simulated halos, expressed as $(1 - f_{b,vir})$, $M_{vir,hyd}$, equates to a mass that is $(1 - f_{b,cosmic})$ times larger than its dmo counterpart. This comparison considers that the dmo counterpart consists solely of a collisionless component, unlike the hydro-simulated halo, which includes both DM and baryonic matter. From Eqs. (1) and (2), it follows that

$$\frac{R_{\Delta,dmo}}{R_{vir,hyd}} = \frac{1 - f_{b,cosmic}}{1 - f_{b,vir}}, \quad (3)$$

and the overdensity $\Delta$ is defined as

$$\Delta = \frac{3 M_{\Delta,dmo}}{4 \pi R_{\Delta,dmo}^3 \rho_c}, \quad (4)$$

with $\rho_c = 3H^3/(8\pi G)$ the cosmic critical density and $H$ the Hubble parameter.

Several authors have studied the validity of the adiabatic approximation for the halo profile (see e.g. Jesseit et al. 2002; Gnedin et al. 2004; Duffy et al. 2010; Velmani & Paranjape 2023), showing that several of the assumptions made are violated; for instance, orbits can deviate strongly from sphericity due to major mergers. Furthermore, energy feedback effects can result in abrupt injection of gas kinetic energy in the cluster core, thus violating the quasi-static assumption. To take these effects into account, we modify Eqs. (2) and (3) as follows:

$$M_{\Delta,dmo} = \frac{1 - f_{b,vir} - \delta_f}{1 - f_{b,cosmic}} M_{vir,hyd}, \quad (5)$$

$$\frac{R_{\Delta,dmo}}{R_{vir,hyd}} = 1 + q \left( \frac{1 - f_{b,cosmic}}{1 - f_{b,vir} - \delta_f} - 1 \right), \quad (6)$$

where $q$ is a quasi-adiabatic parameter that controls the deviation from the adiabatic assumption (see Schneider & Teissier 2015; Paranjape et al. 2021), and $\delta_f$ accounts for the baryonic off-set, for which a halo containing $f_{b,vir} = f_{b,cosmic} - \delta_f$ has the same mass as its dmo counterpart. We note that zero correction corresponds to $\delta_f = 0$ and $q = 1$. It is important to emphasise that these two parameters, $q$ and $\delta_f$, both require calibration against simulations. This calibration is addressed in the following sections.

In summary, to derive the reconstructed DM virial mass, we follow these steps:

1. For a given mass $M_{vir,hyd}$ we calculate the $M_{\Delta,dmo}$ using Eq. (5).

2. We calculate $R_{\Delta,dmo}$ using Eq. (6).

3. The dmo spherical overdensity threshold $\Delta$ is calculated from Eq. (4).

4. The reconstructed dmo virial mass $M_{vir,rec}$ is estimated assuming that the dmo profile follows a Navarro-Frenk-White (NFW, Navarro et al. 1997) profile with a concentration parameter given by Diemer & Joyce (2019).

Therefore, the virial mass $M_{vir,rec}(M_{vir,hyd}, z, f_{b,vir})$ can be reconstructed as a function of the virial mass of the hydro object $M_{vir,hyd}$, its redshift $z$, and its baryonic fraction inside the virial radius $f_{b,vir}$ following the steps outlined above.

## 2.3. Baryonic correction on the halo mass function

The HMF of hydro objects can be expressed as the convolution of the dmo counterpart with the probability distribution function $P(M_{vir,hyd}|M_{vir,mo})$:

$$\frac{dn}{dM_{vir,hyd}} = \int \frac{dn}{dM_{vir,mo}} P(M_{vir,hyd}|M_{vir,mo}) dM_{vir,mo}, \quad (7)$$
where $P(M_{vir,hyd}|M_{vir,dmo})$ is the probability of an object with mass $M_{vir,dmo}$ in the dmo simulation having a counterpart in hydro with mass $M_{vir,hyd}$.

The convolution in Eq. (7) can be simplified if the scatter between $M_{vir,dmo}$ and $M_{vir,hyd}$ is small enough that we can approximate the distribution with a Dirac delta function, as was done in Castro et al. (2020). In this case, we can write

$$ \frac{dn}{dM_{vir,hyd}} = \frac{dn}{dM_{vir,dmo}} \left| \frac{dM_{vir,dmo}}{dM_{vir,hyd}} \right| . \qquad (8) $$

If, in addition, we assume that the baryonic content in halos $f_{b,vir}(M, z)$ of mass $M$ at redshift $z$ is fully predictive, we can make the following substitution in Eq. (8):

$$ M_{vir,dmo} \rightarrow M_{vir,rec}(M_{vir,hyd}, z, f_{b,vir}), \qquad (9) $$

and use the model presented in Sect. 2.2 to compute the HMF in the hydro.

## 2.4. Forecasting Euclid's cluster count observations

As was done by Euclid Collaboration (2023) for the theoretical systematic effects concerning the calibration of the HMF, it is important to quantify how the uncertainties in the baryonic correction model impact the cosmological constraints for a realistic forecast of the cluster sample from the photometric Euclid Wide Survey (Euclid Collaboration 2022). We closely follow the forecasting procedure used in Euclid Collaboration (2021, 2023), and refer to these works for further details. For brevity, here we only highlight the main aspects: synthetic cluster survey data are generated considering a light cone covering 15 000 deg², with redshift range $z \in [0, 2]$ (Laureijs et al. 2011). The number density of halos is sampled assuming the primary calibration for the HMF presented in Euclid Collaboration (2023) with masses corrected by inverting the model described in Sect. 2.2 and assuming that the baryonic content inside the virial overdensity $f_{b,vir}(M_{vir,hyd}, z)$ follows a fiducial relation $f_{b,fid}(M, z)$ extracted from either Magneticum or TNG300 simulations (see Appendix A).

Optical richness $\lambda$ is assigned to the halos of a given virial mass $M_{vir}$ at redshift $z$ according to the relation

$$ \langle \ln \lambda | M_{vir}, z \rangle = \ln A_{\lambda} + B_{\lambda} \ln \left( \frac{M_{vir}}{3 \times 10^{14} h^{-1} M_{\odot}} \right) + C_{\lambda} \ln \left( \frac{E(z)}{E(z=0.6)} \right), \qquad (10) $$

where $E(z)$ is the ratio of the Hubble parameter at redshift $z$ and at redshift 0. We assume a richness range $\lambda \in [20, 2000]$ and a log-normal scatter given by

$$ \sigma^2_{\ln \lambda | M_{vir}, z} = D_{\lambda}^2. \qquad (11) $$

We assume the same fiducial values as in Euclid Collaboration (2023), that is, $A_{\lambda} = 37.8$, $B_{\lambda} = 1.16$, $C_{\lambda} = 0.91$, and $D_{\lambda} = 0.15$, based on the richness-mass relation presented by Saro et al. (2015), and converting it to the virial mass definition assuming that halos follow a NFW profile and the mass-concentration relation given by Diemer & Joyce (2019). The adopted values are in agreement with the results presented by Castignani & Benoist (2016).

Finally, we assume a multivariate Gaussian likelihood with Poisson and sample variance fluctuations given by Hu & Kravtsov (2003), Costanzi et al. (2019),

and Euclid Collaboration (2021). To mimic the uncertainty on the baryon fraction relation and in our baryon correction model, we introduce two nuisance parameters, $\alpha$ and $\beta$. During the cosmological parameter inference, we assume that the baryon fraction is given by

$$ f_{b,vir}(M_{vir,hyd}, z) = (1 + \alpha) f_{b,fid}(M_{vir,hyd}, z), \qquad (12) $$

and, similarly, we assume that the mass reconstruction given by our model relates to the fiducial dmo virial mass as

$$ M_{vir,rec}(M_{vir,hyd}, z, f_{b,vir}) = \frac{M_{vir,dmo}}{1 + \beta}. \qquad (13) $$

We assume Gaussian priors on $\alpha$ and $\beta$ and vary their standard deviation to assess the impact on the cosmological constraints.

# 3. Results

## 3.1. Model calibration

A subset of simulations is leveraged in the calibration of the model, as described in Sect. 2. The choice of simulations is guided by a hierarchical approach aiming to ensure a robust calibration process that incorporates the physical effects incrementally.

For the calibration of the baryonic offset ($\delta_f$), we concentrate on the non-radiative versions of the Magneticum Box 1a. The rationale behind this choice is that these simulations serve as a fundamental calibration point for the baryon fraction, providing a baseline scenario wherein the other effects are assumed to balance to zero.

We used the Magneticum hydrodynamic Box 2 simulations to calibrate the quasi-adiabatic parameter ($q$). This introduces a more advanced level of simulation complexity, building upon the foundational understanding established through the calibration of $\delta_f$.

However, it is important to note that, in principle, both the baryonic offset and the quasi-adiabatic deviation could be influenced by the physics included in the simulations. Therefore, further calibration steps might be required when considering simulations with more detailed subresolution physical models.

### 3.1.1. Baryonic content in non-radiative simulations

In Fig. 1, we present the ratio between the baryonic fraction inside the virial radius and the cosmic baryon fraction, as a function of the virial mass for the three Magneticum simulations for which the non-radiative version is available: C1, C8, and C15 reported in the three columns. The result for the redshifts $z \in \{0.0, 1.0, 2.0\}$ is presented in the different rows. The red lines correspond to the mean and unbiased standard deviations measured for the mass bins. We note that the mean baryonic content inside the virial radius scatters around the following relation despite the background cosmological parameters (marked by the horizontal black line):

$$ \left\langle \frac{f_{b,vir}}{f_{b,cosmic}} \right\rangle = 1 - (0.045 - 0.005 z). \qquad (14) $$

The above relation was obtained by fitting a linear relation for the mean relation as a function of redshift.

Conversely, in Fig. 2, we present the ratio of the virial mass in the non-radiative run with respect to the matched dmo counterpart. Despite the missing baryonic content shown in Fig. 1, we observe that halos have the same average mass in the dmo and
<figure><img src="image_2.png" /><figcaption>Fig. 1. Ratio of the baryonic fraction inside the virial radius with respect to the cosmic baryon fraction as a function of the virial mass for the three Magneticum adiabatic runs C1, C8, and C15 depicted in the three columns. Different rows correspond to the redshifts $z \in \{0.0, 1.0, 2.0\}$. The mean relation, as given by Eq. (14) as a function of redshift, is depicted as the horizontal black line while the red lines correspond to the measured mean and unbiased standard deviations in mass bins. The rightmost column shows the distribution of these ratios for the three simulations.</figcaption></figure>

non-radiative simulations. This leads to the following relation for $\delta_f$ in Eqs. (5) and (6):

$$ \delta_f = (0.045 - 0.005 z) f_{b,\text{cosmic}}. \qquad (15) $$

### 3.1.2. Quasi-adiabatic response

In order to determine the value of the quasi-adiabatic parameter $q$ that appears in Eq. (6) for the model presented in Sect. 2, we proceed in the following way. We bin the halo catalogues of Box 2 of the Magneticum set at five redshifts, $z \in \{0.00, 0.25, 0.50, 1.00, 2.00\}$ in $\log_{10} M_{\text{vir,hyd}}/(h^{-1} M_{\odot})$ with a bin width of 0.3 dex. Within each mass bin, we compute the mean and the unbiased standard deviation of the ratio of the

virial mass in the hydrodynamic and in the dmo simulations, along with the mean of the ratio of the inferred mass and the true dmo mass. Then, the $\chi^2$ between the quasi-adiabatic model and the results from simulations is computed as

$$ \chi^2 = \sum_{i,j} \left( \frac{\langle M_{\text{vir,rec}}/M_{\text{vir,dmo}} \rangle - \langle M_{\text{vir,hyd}}/M_{\text{vir,dmo}} \rangle}{\sigma_{i,j}} \right)_{i,j}^2, \qquad (16) $$

where $i$ and $j$ run over all mass bins and redshifts. In the above equation, $\sigma_{i,j}$ is the error in the binned statistics given by

$$ \sigma_{i,j} = \sqrt{\frac{\text{var}\left(\frac{M_{\text{vir,hyd}}}{M_{\text{vir,dmo}}}\right)_{i,j}}{N_{i,j}} + 0.003^2}, \qquad (17) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Ratio of the virial mass in the non-radiative run with respect to the matched dmo counterpart for the three Magneticum non-radiative runs. Different rows correspond to the redshifts $z \in \{0.0, 1.0, 2.0\}$. The unity relation is depicted as the horizontal black line, while the red lines correspond to the measured mean and unbiased standard deviations in mass bins.</figcaption></figure>

where var is the unbiased variance estimation of the ratio of the virial mass in the hydrodynamic and in the dmo simulations divided by the number of objects in the bin. To this, we add a constant fixed at 0.003 in quadrature. This constant is introduced to adjust for the simplicity of our error modelling, which is estimated directly from the data, possibly understating the true uncertainties in the relation. The addition of this constant ensures a resulting $\chi^2$ of the order of the unity.

In Fig. 3, we plot the values of the $\Delta\chi^2$ as a function of $q$, both for single redshifts and for the analysis based on combining all redshifts. The $\Delta\chi^2$ is defined as:

$$ \Delta\chi^2 = \chi^2 - \chi_{\text{best-fit}}^2, \qquad (18) $$

where $\chi_{\text{best-fit}}^2$ corresponds to the best fit of the joint analysis. We observe that all redshifts prefer significant deviations from the adiabatic prediction ($q = 1$). All redshifts present a minimum around the joint analysis minimum $q = 0.373$ (presented as the

vertical black dotted line). Exceptions occur at $z = 0.25$ and $z = 2.0$ that have their best fit shifted to $q \approx 0.23$ and $q \approx 0.60$, respectively. Still, assuming Gaussianity, those shifts correspond to only a $2\sigma$ deviation for $z = 0.25$ and less than $1\sigma$ for $z = 2.0$. Given the lack of strong statistical significance, we interpret this deviation as likely a result of statistical variation rather than a reflection of any substantive physical process at this particular redshift. Therefore, in the following, we assume the joint analysis best fit for $q$ to hold at all redshifts.

## 3.2. Validation of the model

In Fig. 4, we present the ratios of the original halo masses in hydro simulations (top panel) and the reconstructed halo masses (bottom panel), both to the dmo counterpart. Results are shown for the subset of Magneticum Boxes 1a that assume different cosmological parameters (see Table 2), at four redshifts, $z \in$
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Variation of the $\chi^2$ with respect to its minimum, as a function of the quasi-adiabatic parameter $q$ for the model presented in Sect. 2. We present the results for the joint analysis of five redshifts $z \in \{0.00, 0.25, 0.50, 1.00, 2.00\}$ from the Magneticum Box 2 simulation, as well as the results at each redshift. The best fit for $q = 0.373$ is marked with the vertical dashed line.</figcaption></figure>

$\{0.0, 0.3, 0.5, 0.9\}$. The medians of the corresponding distributions are shown in Fig. 5. The bottom panel presents the model's performance in recovering the dmo mass for the simulations with different cosmologies colour-coded by the baryon fraction.

As expected, the median value of the ratio between hydro and dmo halo masses shows significant deviations from unity, in a way that becomes more significant as the value of the baryon fraction increases. On the other hand, we see from the bottom panels that our model, on average, allows us to correctly recover the halo masses from the dmo simulations. Quite remarkably, this accuracy is independent of the cosmological model adopted in the simulations, as the correlation of the net effect with the baryonic fraction is absent in the reconstructed mass. In general, dmo masses are recovered by our model with subpercent accuracy. Possible exceptions are represented by the simulations with the largest cosmic baryon fraction at $z = 0$, for which our model seems to slightly over-predict the reconstructed mass by roughly 1 percent. Lastly, we observe that subdominant to the model target accuracy of 1 percent, there is a correlation between the model performance and the redshift, with the median of the ratio decreasing with increasing redshift. This suggests that there is a feature missing in our model, such as a redshift-dependent quasi-adiabatic factor reflecting that at higher redshift, the effect is driven by radiative cooling, while at low redshift AGN feedback is the main agent (see Castro et al. 2020). An increase in the sophistication of the model to push our model accuracy to the subpercent level is left for further investigation in future work.

Similarly to Fig. 4, in Fig. 6 we present the performance of the model for the simulations with C8 (WMAP7) background cosmology, but assuming different values for the parameters that define the efficiency of feedback from both AGN and SN. Vertical lines denote the median of the corresponding distribution. In particular, we varied (a) the efficiency of the AGN feedback by changing the value of the BlackholeFeedbackFactor parameter from 0.1 to 0.2 (default value was 0.15), and (b) the value of the velocity of the SN-driven galactic outflows in the model by Springel & Hernquist (2003) from 350 km s⁻¹ to 500 km s⁻¹ and 800 km s⁻¹, respectively. In general, we note that the effect of changing the parameters regulating stellar and AGN feedback is relatively small, especially at low redshift. In any case, the accuracy with which our model recovers the halo masses in the dmo

simulation is not degraded by changing the feedback parameters. The difference between the vertical lines is always smaller for the reconstructed mass than when considering the actual masses from the hydro simulations.

Similarly to Figs. 4 and 6, in Figs. 7 and 8 we present the performance of the model on the Magneticum Box 2b and on the TNG300 simulations, respectively. The different columns present the results for $z \in \{0.2, 0.5, 1.0, 2.0\}$ for Magneticum and $z \in \{0.0, 0.5, 1.0, 2.0\}$ for TNG300. We note that considering Box 2b instead of Box 1a of the Magneticum set allows us to validate the robustness of our method when increasing resolution. At the same time, the analysis of the TNG300 box allows us to further stress the validity of our model at an even higher resolution and for a substantially different hydrodynamical solver, star-formation model, and implementation of both stellar and AGN feedback. We present the results for two mass regimes: $3 \times 10^{13} \le M_{\text{vir,dmo}}/(h^{-1}M_{\odot}) < 10^{14}$ and $M_{\text{vir,dmo}}/(h^{-1}M_{\odot}) \ge 10^{14}$. We refer to the former as the group regime and the latter as the regime of galaxy clusters. We note that, due to the smaller size of the TNG300 box, this simulation does not contain cluster-sized halos at $z = 2$. Both simulations predict a stronger impact on the group regime than on the cluster regime; this is expected as the potential well is shallower in the former regime than in the latter, and therefore the baryonic depletion requires less feedback. In general, TNG300 predicts a smaller impact of baryonic effects on halo masses in the cluster regime than the Magneticum Box 2b simulation. In any case, our model performs equally well in both simulations and mass regimes. The smaller impact predicted by the TNG300 simulations is due to the fact that Magneticum and TNG300 simulations predict a similar baryonic content for virial halos at $z = 1$ as can be seen in Fig. A.1. However, the cosmic baryon fraction assumed in the Magneticum reference simulation is ~7 percent higher, and therefore it can be inferred that the Magneticum simulations exhibit a more pronounced gas depletion and more active feedback mechanisms compared to those in TNG. This results in a greater impact on halo mass in Magneticum due to a more significant displacement of gas.

### 3.3. Robustness of the model to the assumed baryon fraction relation

Above, we demonstrate the performance of our model in reconstructing $M_{\text{vir,dmo}}$ under the assumption that the baryon fraction for each halo is known. Additionally, it is of interest to evaluate the performance of the model when employing an averaged baryon fraction-mass relation across a sample, as an alternative approach. This approach is arguably more realistic for implementation in cosmological analyses.

In Fig. 9, we present the ratio between the cumulative halo abundance assuming different masses in the hydro simulation and the halo abundance of virial halos in the dmo simulation. We present the results for the virial mass of the hydro simulation as well as the virial prediction of our model when taking into account the following assumptions:

1. That the baryon fraction of individual objects is known.

2. A predictive baryon fraction relation given by Eq. (A.1) with parameters shown in Table A.1 for the Magneticum set.

3. A stochastic baryon fraction with mean given by Eq. (A.1) with parameters shown in Table A.1 for the Magneticum set and with a scatter of 0.17 dex.

The third case addresses the possibility that the scatter of the cluster gas fraction in the simulation is underestimated compared to observational data. We assume the scatter of
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Ratio between either the hydrodynamic mass (top panels) or the reconstructed virial mass (bottom panels) and the dmo counterpart on the subset of Magneticum Boxes 1a for the simulations with different cosmologies (see Table 2). Different columns correspond to the four redshifts $z \in \{0.0, 0.3, 0.5, 0.9\}$ and different colours correspond to different cosmic baryon fractions $f_b$. The filled dark-grey and light-grey areas mark the 1 and 2 percent halo mass variation induced by baryonic effects, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Median of the distributions shown in Fig. 4 as a function of the cosmic baryon fraction $f_b$ for the ratio between either the hydrodynamic mass (dotted lines) or the reconstructed virial mass (full lines) and the dmo counterpart on the subset of Magneticum Boxes 1a for the simulations with different cosmologies (see Table 2). Different colours correspond to the four redshifts $z \in \{0.0, 0.3, 0.5, 0.9\}$. The filled regions correspond to one median absolute deviation around the median.</figcaption></figure>

0.17 dex reported by Andreon et al. (2017) for the gas fraction. To increase the readability of the plot, we only present error bars for the first case, assuming uncorrelated Poisson errors for the abundances. Different panels correspond to the redshifts $z \in \{0.0, 0.5, 1.0, 2.0\}$, and for all panels we use Magneticum Box 2b, except for $z = 0.0$, where we use Box 2, as Box 2b was not run to this redshift.

In Fig. 9, the baryonic physics tends to produce halos that are less massive than their dmo counterparts, as previously discussed by Castro et al. (2020). Therefore, at fixed mass value,

the abundance of halos in a hydro simulation is lower than that derived from a dmo simulation (blue lines). Conversely, converting $M_{vir,hyd}$ to the corresponding dmo virial masses, the abundances of the two simulations match (green and orange lines). Regarding the performance of our model, we note that using a collective relation (red and green lines) instead of the individual object baryon fraction (orange line) does not significantly degrade the mass reconstruction. This is neither obvious nor expected a priori, as non-linear functions do not commute with the median operation. The robustness of the performance comes from the tightness of the baryon fraction relation at fixed mass (see Fig. A.1), and ensures the applicability of our model in cluster cosmology studies, even if we assume a larger variance for the baryon fraction distribution. Therefore, the approximations that lead to the simplification of Eq. (7)–(8) are validated and do not statistically affect the model’s performance.

**4. Impact of baryonic effects on a Euclid-like cluster abundance analysis**

In this section, we quantify the biases in deriving the posteriors on the cosmological parameters $\Omega_m$ and $\sigma_8$ caused by neglecting the baryonic effects on halo masses. Specifically, we assess these biases using a fully predictive baryon fraction relation measured from either the Magneticum or the TNG300 simulations when creating the synthetic catalogue (see Appendix A for more details). Following this, we evaluate the efficacy of our method in recovering halo masses in dmo simulations.

We perform the forecast following the methodology described in Sect. 2.4. For the likelihood analysis, we assume flat priors on the cosmological parameters and Gaussian priors on the mass-observable parameters of Eqs. (10)
<figure><img src="image_7.png" /><figcaption>Fig. 6. Similar to Fig. 4 but for the simulations with C8 (WMAP7) background cosmology that assume different values of the parameters describing SN and AGN feedback.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Ratio of the hydrodynamic mass (top panels) and of the reconstructed virial mass (bottom panels) to the dmo counterpart on the Magneticum Boxes 2b. The different columns present the results for $z \in \{0.2, 0.5, 1.0, 2.0\}$. We show the results for two mass regimes: $3 \times 10^{13} \le M_{vir,dmo}/(h^{-1} M_{\odot}) < 10^{14}$ (blue histograms) and $M_{vir,dmo}/(h^{-1} M_{\odot}) \ge 10^{14}$ (orange histograms).</figcaption></figure>

and (11) with the mean given by the fiducial values and with a rms of 1 and 3 percent, respectively. The likelihood sampling is performed following the ensemble slice sampling method (see Karamanis & Beutler 2021) implemented in ZEUS (see Karamanis et al. 2021).

Firstly, we assess the tension between the inferred cosmological parameters and the fiducial ones after marginalising over

the other parameters, as quantified by the index of inconsistency (IOI; Lin & Ishak 2017), which is calculated as

$$IOI = \frac{\delta^t \Sigma^{-1} \delta}{2} \qquad (19)$$

In the above expression, $\delta$ is the two-dimensional difference vector between the best-fit values and the fiducial values of the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Similar to Fig. 7 but for the TNG300 simulation.</figcaption></figure>

$\Omega_m$ and $\sigma_8$ cosmological parameters, while $\Sigma$ is the covariance matrix between these parameters, which we assume to be Gaussian distributed. We calculate the IOI when ignoring the baryonic impact, assuming the cosmological parameters and the virial baryon content relation from the TNG and the Magneticum simulations.

In Table 3, we present the summary statistics for the forecast of the impact of the treatment of the baryonic effects on halo masses. We present the results for the constraints on $\Omega_m$ and $\sigma_8$ to be obtained from the cluster counts analysis of the forthcoming Euclid survey. Ideally, the IOI should be kept below 1 to ensure that any correlations with other systematic effects do not amplify tensions in the final results (Euclid Collaboration 2024).

For both Magneticum and TNG simulations, we observe an increasing impact of the baryonic effects on halo masses as stronger priors on the mass-observable scaling relation are assumed. This is an expected result, as a tighter prior on the scaling relation results in tighter constraints on the cosmological parameters, making cosmological posteriors more sensitive to systematic effects. While ignoring baryonic effects in the Magneticum simulations would lead to an IOI of 4.1 and 2.0 for 1 and 3 percent priors, respectively, for TNG300 these values reduce to 1.8 and 0.9.

The reasons for the smaller impact of baryons predicted by TNG with respect to Magneticum are two fold.

1. Baryon content in TNG300: As mentioned in Sect. 3.2, the baryon content in massive clusters observed in TNG300 more closely aligns with the cosmic baryon fraction assumed in this simulation. This alignment suggests a smaller depletion of baryons due to feedback mechanisms in TNG300 compared to Magneticum.

2. Behaviour of scaling relation parameters: In the Magneticum simulations, the recovered posteriors for the scaling relation parameters closely match the fiducial values. However, in the TNG simulations, we observe a notable shift toward lower values – of approximately $1\sigma$ – in the posteriors for $C_\lambda$, that is, the parameter that controls redshift evolution.

This shift indicates that the scaling relation in TNG absorbs some of the baryonic impacts. Such absorption by $C_\lambda$ is somewhat expected because it is the sole parameter in our mass-observation relation that evolves with background evolution, similar to the baryon fraction. The specific redshift dependency of the baryon fraction as a function of mass for TNG300 is detailed in Table A.1.

We note that, in the Bocquet et al. (2019) cluster cosmology analysis from the SPT-SZ Survey, $C_\lambda$ is especially degenerate with the dark energy equation of state parameter $w$, which we have not considered as a cosmological parameter for our forecast. The absorption of the baryonic impact on $C_\lambda$ is likely to lead to biased dark energy constraints in this case.

In Table 3, we also report the degradation in the constraints on cosmological parameters after also marginalizing over the limited knowledge of the baryon content and the model accuracy, which are described by the $\alpha$ and $\beta$ parameters (see Eqs. (12) and (13)). We present the relative change in the figure of merit (FOM; see Huterer & Turner 2001; Albrecht et al. 2006) in the ($\Omega_m, \sigma_8$) plane for different scenarios concerning hypothetical perfect knowledge of the baryonic content and a perfect performance of the baryonic correction model. Again, we consider 1 and 3 percent priors on the richness-mass scaling relations. In all cases, we assume a Gaussian prior on $\beta$ with zero mean and 0.01 standard deviation in agreement with the accuracy and precision of our model. Concerning the priors on $\alpha$, we consider three scenarios: 0 (perfect knowledge), 1 percent, and 10 percent. To provide some context, Chiu et al. (2018) claim that the characteristic baryon mass fraction inside $R_{500c}$ is about $(12.8 \pm 1.3)$ percent for clusters with a mass of $M_{500c} = 4.8 \times 10^{14} M_\odot$ and at redshift $z = 0.6$ based on their analysis of a set of 91 galaxy clusters from the 2500 deg$^2$ South Pole Telescope SPT-SZ survey (Bleem et al. 2015). Therefore, our value of 10 percent for the prior on $\alpha$ is in line with the current level of precision in the observational calibration of this relation. For 1 and 3 percent priors in the scaling relation and perfect knowledge of the baryon fraction relation, we observe that marginalizing
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Ratio of the cumulative halo abundance in the hydro simulation to the halo abundance of virial halos in the dmo simulation. The results are shown for the simulated virial mass and the prediction of our model based on the individual baryon fractions, a predictive baryon fraction relation, and a stochastic baryon fraction relation. Error bars are displayed for the first case, assuming uncorrelated Poisson errors. Different panels correspond to redshifts $z \in \{0.0, 0.5, 1.0, 2.0\}$ using Magneticum Box 2b, except for $z = 0.0$, where Box 2 is used. The grey regions correspond to 2.5 and 5 percent.</figcaption></figure>

over the accuracy of our model impacts the FOM in the cosmological constraints by 15 and 11 percent, respectively. While 1 percent prior in the baryon fraction relation does not significantly further degrade the FOM for the case where the richness-mass relation is known at 3 percent, it reduces the FOM by a further

7 percent if the richness-mass relation is known with 1 percent precision. Lastly, for a 10 percent prior in the baryon fraction relation, the impact in the FOM is roughly 70 and 40 percent for the cases where the richness-mass relation is known at 1 and 3 percent precision level, respectively.

For all cases considered for the forecast, the correlation between $\alpha$ and $\beta$ and the cosmological parameters is less significant than the correlation between the cosmological parameters themselves. Although the absence of strong correlations between cosmological parameters and $\alpha$ and $\beta$ is a desired feature in order to improve the robustness of the analysis, in the specific case of $\alpha$, it is also due to the simplified functional form considered for the uncertainty on the baryon fraction relation. A more detailed analysis of the impact of the baryon fraction relation, considering a more flexible parametrisation, is out of the scope of this paper and is left for further investigation.

In summary, the results presented in this section confirm that neglecting baryonic effects on halo masses leads to significantly biased constraints on cosmological parameters from the cluster number counts, especially when strong priors on the mass-observable relation are assumed. On the other hand, we also demonstrate that our model, which is designed to correct for such baryonic effects, significantly reduces this bias.

## 5. Conclusions

In this paper, we present a model of the impact of baryonic feedback on the virial masses of galaxy clusters. The main aim of our analysis is to verify whether such baryonic effects on halo masses can be modelled with sufficient precision that they do not represent a limiting factor in the cosmological constraining power of the survey of galaxy clusters that will be obtained during the *Euclid* wide photometric survey.

Our model assumes that the feedback effects on the intracluster gas show a quasi-adiabatic behaviour and have the baryon fraction inside the virial radius as input. In this case, this effect can be reliably described by a phenomenological model that depends only on two free parameters and can be calibrated against cosmological hydro simulations. These two parameters control the minimum baryonic depletion observed in non-radiative simulations and the deviation from the adiabatic behaviour. Although our model is calibrated using a set of non-radiative hydrodynamic simulations and a single realisation of a full-physics simulation, we demonstrate that its performance is resilient to changes in the subresolution physics and cosmological parameters. This result is established by confronting our model with different simulations from the Magneticum suite and the TNG300 simulation. Finally, with the resulting baryonic correction model, we assess the impact of its accuracy and precision in the cosmological parameters inference from an idealised *Euclid* cluster number count experiment. Our main conclusions can be summarised as follows.

- We show that the baryonic feedback effects on the intracluster gas can be accurately modelled using a quasi-adiabatic approach with two parameters: one controlling the minimum baryonic depletion and the other the deviation from adiabatic behaviour.

- We successfully calibrated these two parameters against cosmological hydro simulations. This model calibration was conducted using a set of non-radiative hydrodynamic simulations as well as a single full-physics simulation.

- Our model attains one percent relative accuracy in determining the virial dmo equivalent mass of simulated galaxy clusters.
<table><caption>**Table 3.** Summary statistics for the forecast of the impact of the baryon mass correction model on Euclid's cluster count cosmological constraints on $\Omega_m$ and $\sigma_8$.</caption><thead><tr><th rowspan="2">Summary statistics</th><th rowspan="2">Richness-mass relation priors</th><th rowspan="2">Baryonic treatment</th><th colspan="2">⟨f<sub>b,vir</sub>(M, z)⟩</th><th rowspan="2">Value</th></tr><tr><th>Relation</th><th>Prior on α</th></tr></thead><tbody><tr><td rowspan="2">IOI</td><td rowspan="2">1%</td><td rowspan="2">Ignoring</td><td>TNG300</td><td>–</td><td>1.8 ± 0.1</td></tr><tr><td>Magneticum</td><td>–</td><td>4.1 ± 0.2</td></tr><tr><td rowspan="2">FOM<sub>corrected</sub>/FOM</td><td rowspan="2">3%</td><td rowspan="2">Ignoring</td><td>TNG300</td><td>–</td><td>0.9 ± 0.1</td></tr><tr><td>Magneticum</td><td>–</td><td>2.0 ± 0.2</td></tr><tr><td rowspan="4">FOM<sub>corrected</sub>/FOM</td><td rowspan="2">1%</td><td rowspan="2">Our model</td><td rowspan="2">Magneticum</td><td>Perfect knowledge</td><td>0.85 ± 0.02</td></tr><tr><td>1%</td><td>0.78 ± 0.02</td></tr><tr><td rowspan="2">3%</td><td rowspan="2">Our model</td><td rowspan="2">Magneticum</td><td>Perfect knowledge</td><td>0.89 ± 0.03</td></tr><tr><td>1%</td><td>0.90 ± 0.03</td></tr><tr><td></td><td></td><td></td><td></td><td>10%</td><td>0.61 ± 0.02</td></tr></tbody></table>

**Notes.** The IOI quantifies the tension in the posteriors for the different treatments with respect to hypothetical perfect knowledge of the model used to create the synthetic data (see text). The FOM ratio assesses the attenuation of the constraining power of the cluster counts once one marginalises the uncertainties of the model and the baryon content of the clusters. For both statistics, errors were estimated by sampling multivariate Gaussian distributions with the correlation extracted from the respective MCMC chains and recomputing the statistics from them.

- The robustness of our model is demonstrated using different simulations from the Magneticum suite and the TNG300 simulation. Our findings indicate that the model's effectiveness does not significantly deviate with changes to the subresolution physics or cosmological parameters, showing strong consistency between our model predictions and these independent simulations.

- Unlike previous works (Cui et al. 2014; Bocquet et al. 2016; Castro et al. 2020), our model does not rely on a single baryonic physics model or assume the universality of the HMF as a function of cosmology. This flexibility represents a significant advantage over other calibrations of the impact of baryonic feedback on the HMF, which cannot predict outcomes under altered baryonic models and implicitly assume HMF universality.

- The findings of our research substantiate previous claims concerning the potentially significant impact of baryons – if neglected – on the cosmological constraints derived from the Euclid photometric galaxy cluster number counts.

- Most importantly, the uncertainties linked to our model for correcting baryonic effects on cluster masses are shown to be subdominant to both the precision of the expected calibration of the mass-observable relation with Euclid and our current understanding of the baryon fraction within galaxy clusters. Lastly, it is important to note that the model presented in this paper is based on the assumption that AGN feedback acts in a quasi-adiabatic manner. Future work could extend the current framework by investigating its validity concerning other prescriptions for AGN feedback as well as additional non-thermal processes not explored in this suite of simulations.

**Acknowledgements.** TC, SB, AS and AF are supported by the INFN INDARK PD51 grant. TC, SB and AS are supported by the Fondazione ICSC National Recovery and Resilience Plan (PNRR) Project ID CN-00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: “Potenziamento strutture di ricerca e creazione di “campioni nazionali” di R&S (M4C2-19)” – Next Generation EU (NGEU). TC and AS are also supported by the FARE MIUR grant ‘ClustersXEuclid’ R165SBKTMA. AS is also supported by the ERC ‘ClustersXCosmo’ grant agreement 716762. KD acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311 as well as support through the COMPLEX project from the European Research Council

(ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. AR is supported by the PRIN-MIUR 2017 WSCC32 ZOOMING grant. We acknowledge the computing centre of CINECA and INAF, under the coordination of the “Accordo Quadro (MoU) per lo svolgimento di attività congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione”, for the availability of computing resources and support. AF acknowledges support from Brookhaven National Laboratory. We acknowledge the use of the HOTCAT computing infrastructure of the Astronomical Observatory of Trieste – National Institute for Astrophysics (INAF, Italy) (see, Bertocco et al. 2020; Taffoni et al. 2020). The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d'Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org).

## References

- Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, PRD, 102, 023509

- Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv e-prints [arXiv:astro-ph/0609591]

- Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409

- Andreon, S., Wang, J., Trinchieri, G., Moretti, A., & Serra, A. L. 2017, A&A, 606, A24

- Aricò, G., Angulo, R. E., Contreras, S., et al. 2021, MNRAS, 506, 4070

- Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110

- Bertocco, S., Goz, D., Tornatore, L., et al. 2020, ASP Conf. Ser., 527, 303

- Bleem, L. E., Stalder, B., De Haan, T., et al. 2015, ApJS, 216, 27

- Blumenthal, G. R., Faber, S. M., Flores, R., et al. 1986, ApJ, 301, 27

- Bocquet, S., Saro, A., Mohr, J. J., et al. 2015, ApJ, 799, 214

- Bocquet, S., Saro, A., Dolag, K., et al. 2016, MNRAS, 456, 2361

- Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55

- Bocquet, S., Heitmann, K., Habib, S., et al. 2020, ApJ, 901, 5

- Bondi, H. 1952, MNRAS, 112, 195

- Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273

- Borgani, S., & Kravtsov, A. 2011, Adv. Sci. Lett., 4, 204

- Borgani, S., Rosati, P., Tozzi, P., et al. 2001, ApJ, 561, 13

- Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80

- Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343

- Castignani, G., & Benoist, C. 2016, A&A, 595, A111

- Castro, T., Quartin, M., Giocoli, C., et al. 2018, MNRAS, 478, 1305

- Castro, T., Borgani, S., Dolag, K., et al. 2020, MNRAS, 500, 2316
Chiu, I., Mohr, J. J., McDonald, M., et al. 2018, MNRAS, 478, 3072

Costanzi, M., Rozo, E., Simet, M., et al. 2019, MNRAS, 488, 4779

Costanzi, M., Saro, A., Bocquet, S., et al. 2021, PRD, 103, 043522

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Cui, W., Borgani, S., & Murante, G. 2014, MNRAS, 441, 1769

Davis, M., Efstathiou, G., Frenk, C. S., et al. 1985, ApJ, 292, 371

Debackere, S. N. B., Schaye, J., & Hoekstra, H. 2021, MNRAS, 505, 593

Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486

Di Matteo, T., Colberg, J., Springel, V., et al. 2008, ApJ, 676, 33

Diemer, B., & Joyce, M. 2019, ApJ, 871, 168

Dolag, K., Borgani, S., Murante, G., et al. 2009, MNRAS, 399, 497

Duffy, A. R., Schaye, J., Kay, S. T., et al. 2010, MNRAS, 405, 2161

Ellien, A., Durret, F., Adami, C., et al. 2019, A&A, 628, A34

Euclid Collaboration (Adam, R., et al.) 2019, A&A, 627, A23

Euclid Collaboration (Fumagalli, A., et al.) 2021, A&A, 652, A21

Euclid Collaboration (Scaramella, R., et al.) 2022, A&A, 662, A112

Euclid Collaboration (Castro, T., et al.) 2023, A&A, 671, A100

Euclid Collaboration (Deshpande, A. C., et al.) 2024, A&A, 684, A138

Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761

Fumagalli, A., Costanzi, M., Saro, A., et al. 2023, A&A, 682, A148

Gnedin, O. Y., Kravtsov, A. V., Klypin, A. A., et al. 2004, ApJ, 616, 16

Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, JCAP, 07, 008

Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304

Holder, G., Haiman, Z., & Mohr, J. 2001, ApJ, 560, L111

Hoyle, F., & Lyttleton, R. A. 1939, Math. Proc. Camb. Philos. Soc., 35, 405

Hu, W., & Kravtsov, A. V. 2003, ApJ, 584, 702

Huterer, D., & Turner, M. S. 2001, PRD, 64, 123527

Jesseit, R., Naab, T., & Burkert, A. 2002, ApJ, 571, L89

Karamanis, M., & Beutler, F. 2021, Stat. Comput., 31, 61

Karamanis, M., Beutler, F., & Peacock, J. A. 2021, MNRAS, 508, 3589

Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2011, ApJS, 192, 18

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

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints
[arXiv:1110.3193]

Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022, A&A, 659, A88

Lin, W., & Ishak, M. 2017, PRD, 96, 023532

Mantz, A. B., Von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205

Martizzi, D., Teyssier, R., Moore, B., et al. 2013, MNRAS, 432, 1947

McCarthy, I. G., Schaye, J., Bird, S., et al. 2017, MNRAS, 465, 2936

McDonald, M., Bayliss, M., Benson, A., et al. 2012, Nature, 488, 349

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

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

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2

Pakmor, R., & Springel, V. 2013, MNRAS, 432, 176

Pakmor, R., Bauer, A., & Springel, V. 2011, MNRAS, 418, 1392

Pakmor, R., Marinacci, F., & Springel, V. 2014, ApJ, 783, L20

Paranjape, A., Choudhury, T. R., & Sheth, R. K. 2021, MNRAS, 503, 4147

Pillepich, A., Springel, V., Nelson, D., et al. 2018a, MNRAS, 473, 4077

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018b, MNRAS, 475, 648

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

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

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

Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645

Ryden, B. S., & Gunn, J. E. 1987, ApJ, 318, 15

Saro, A., Bocquet, S., Rozo, E., et al. 2015, MNRAS, 454, 2305

Sartoris, B., Biviano, A., Fedeli, C., et al. 2016, MNRAS, 459, 1764

Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schaye, J., Kugel, R., Schaller, M., et al. 2023, MNRAS, 526, 4978

Schellenberger, G., David, L., O'Sullivan, E., Vrtilek, J. M., & Haines, C. P. 2019, ApJ, 882, 59

Schneider, A., & Teyssier, R. 2015, JCAP, 12, 049

Singh, P., Saro, A., Costanzi, M., et al. 2020, MNRAS, 494, 3728

Springel, V. 2005, MNRAS, 364, 1105

Springel, V. 2010, MNRAS, 401, 791

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Springel, V., White, S. D. M., Tormen, G., et al. 2001a, MNRAS, 328, 726

Springel, V., Yoshida, N., & White, S. D. M. 2001b, New Astron., 6, 79

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Springel, V., Pakmor, R., Zier, O., et al. 2021, MNRAS, 506, 2871

Steigman, G., Sarazin, C. L., Quintana, H., et al. 1978, AJ, 83, 1050

Taffoni, G., Becciani, U., Garilli, B., et al. 2020, ASP Conf. Ser., 527, 307

Teyssier, R., Moore, B., Martizzi, D., et al. 2011, MNRAS, 414, 195

Tinker, J. L., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Tornatore, L., Borgani, S., Dolag, K., et al. 2007, MNRAS, 382, 1050

van Daalen, M. P., Schaye, J., Booth, C. M., & Vecchia, C. D. 2011, MNRAS, 415, 3649

van Daalen, M. P., McCarthy, I. G., & Schaye, J. 2020, MNRAS, 491, 2424

Velliscig, M., van Daalen, M. P., Schaye, J., et al. 2014, MNRAS, 442, 2641

Velmani, P., & Paranjape, A. 2023, MNRAS, 520, 2867

Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Vogelsberger, M., Marinacci, F., Torrey, P., et al. 2020, Nat. Rev. Phys., 2, 42

Watson, W. A., Iliev, I. T., D'Aloisio, A., et al. 2013, MNRAS, 433, 1230

Webb, T., Noble, A., DeGroot, A., et al. 2015, ApJ, 809, 173

Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291

Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99

Yuan, T., Elagali, A., Labbé, I., et al. 2020, Nat. Astron., 4, 957

Zeldovich, Y. B., Klypin, A., Khlopov, M. Y., et al. 1980, Sov. J. Nucl. Phys., 31, 664

1 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy
e-mail: tiago.batalha@inaf.it

2 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, TS, Italy

3 IFPU, Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste, Italy

4 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy

5 Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy

6 Institute for Computational Science, University of Zurich, Winterhuterstrasse 190, 8057 Zurich, Switzerland

7 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany

8 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy

9 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Via Piero Gobetti 93/2, 40129 Bologna, Italy

10 Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France

11 Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France

12 Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK

13 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy

14 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

15 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

16 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany

17 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, (TO), Italy

18 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy

19 INFN-Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

20 Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126 Napoli, Italy

21 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

22 INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy

23 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

24 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy

25 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

26 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy

27 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, (Barcelona), Spain
28 Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra, (Barcelona), Spain

29 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany

30 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy

31 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

32 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

33 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK

34 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

35 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain

36 University of Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, 69622 Villeurbanne, France

37 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

38 UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France

39 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

40 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande 1749-016, Lisboa, Portugal

41 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande 1749-016, Lisboa, Portugal

42 Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

43 INAF-Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere, 100, 00100 Roma, Italy

44 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK

45 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

46 INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

47 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany

48 Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern 0315, Oslo, Norway

49 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

50 von Hoerner & Sulger GmbH, SchloßPlatz 8, 68723 Schwetzingen, Germany

51 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark

52 Cosmic Dawn Center (DAWN), Copenhagen, Denmark

53 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

54 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 Helsinki, Finland

55 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

56 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, 91191 Gif-sur-Yvette, France

57 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland

58 Department of Physics, PO Box 64, 00014 University of Helsinki, Finland

59 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

60 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands

61 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

62 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

63 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK

64 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France

65 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

66 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France

67 Space Science Data Center, Italian Space Agency, Via del Politecnico snc, 00133 Roma, Italy

68 Centre National d'Etudes Spatiales - Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

69 Institute of Space Science, Str. Atomistilor, nr. 409 Măgurele, Ilfov 077125, Romania

70 Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, 38204 San Cristóbal de La Laguna, Tenerife, Spain

71 Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain

72 Dipartimento di Fisica e Astronomia "G. Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy

73 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

74 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile

75 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria

76 Institut d'Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2-4, 08034 Barcelona, Spain

77 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain

78 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain

79 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

80 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal

81 Universidad Politécnica de Cartagena, Departamento de Elétrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain

82 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France

83 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

84 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy

85 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

86 Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France

87 CEA Saclay, DFR/IRFU, Service d'Astrophysique, Bat. 709, 91191 Gif-sur-Yvette, France

88 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany

89 Université St Joseph, Faculty of Sciences, Beirut, Lebanon

90 Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France

91 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France

92 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, TS, Italy

93 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain

94 CERCA/ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA
95 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

96 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

97 NASA Ames Research Center, Moffett Field, CA 94035, USA

98 Kavli Institute for Particle Astrophysics & Cosmology (KIPAC),
Stanford University, Stanford, CA 94305, USA

99 Bay Area Environmental Research Institute, Moffett Field, CA 94035, USA

100 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA

101 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy

102 Institute Lorentz, Leiden University, PO Box 9506, Leiden 2300 RA, The Netherlands

103 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

104 Department of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA

105 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain

106 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada

107 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

108 Department of Computer Science, Aalto University, PO Box 15400, Espoo 00 076, Finland

109 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany

110 Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

111 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000 Grenoble, France

112 Department of Physics and Astronomy, Veskilinnantie 5, 20014 University of Turku, Finland

113 Serco for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692 Madrid, Spain

114 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Victoria 3122, Australia

115 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia

116 Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm 106 91, Sweden

117 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK

118 Centre de Calcul de l'IN2P3/CNRS, 21 avenue Pierre de Coubertin, 69627 Villeurbanne Cedex, France

119 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy

120 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 – c/o Dipartimento di Fisica, Edificio G. Marconi 00185 Roma, Italy

121 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal

122 Zentrum für Astronomie, Universität Heidelberg, Philosophenweg 12, 69120 Heidelberg, Germany

123 Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, Roma, Italy

124 INFN, Sezione di Roma 2, Via della Ricerca Scientifica 1, Roma, Italy

125 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333, CA Leiden, The Netherlands

126 Department of Physics and Astronomy, University of California, Davis, CA 95616, USA

127 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA

128 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark
## Appendix A: Baryon content of galaxy clusters in full-physics simulations

In this Appendix, we report results pertaining to the baryon content of galaxy clusters and groups within the halo virial radius for both the Magneticum and the TNG300 simulations. The results provided here are used in our analysis for the calibration of the parameters that define our model to correct halo masses for baryonic effects. A comparison with observational data is beyond the scope of the present study.

In Fig. A.1 we show the baryonic fraction inside the virial halos as a function of the halo virial mass for both the Magneticum (upper panels) and TNG300 (lower panels) simulations, at four different redshifts $z \in \{0.0, 0.5, 1.0, 2.0\}$. The cosmic baryon fraction assumed by the respective simulations is shown as the horizontal black line. The median baryon fraction in halos is shown by the red line and the best fit for the relation is shown in green. For the Magneticum results, we show the results obtained for Box 2 and Box 1a with orange and blue symbols, respectively. Box 2 samples the mass range of groups and low-mass clusters at increased resolution, while the larger size of Box 1a allows us to sample the most massive clusters at lower resolution (see Table 1).

The best-fit relations for the Magneticum and TNG300 are given by

<table><caption>**Table A.1.** Best-fit parameters for the fitting function presented in Eq. (A.1) for the baryon fraction inside the virial radius for the Magneticum and TNG300 simulations.</caption><thead><tr><th>Simulation</th><th>γ(z)</th><th>δ(z)</th></tr></thead><tbody><tr><td>TNG300</td><td>0.4106 z − 1.4640</td><td>0.0397 z<sup>2</sup> − 0.04358 z + 0.03567</td></tr><tr><td>Magneticum</td><td>0.7008 z − 1.7505</td><td>0.2</td></tr></tbody></table>

$$f_{\text{b,fid}}(m, z) = f_{\text{b,cosmic}} m^{-\gamma} (1 + m^{1/\delta})^{\gamma\delta}, \qquad (\text{A.1})$$

with m given by

$$m = \log_{10} \frac{M_{\text{vir}}}{10^{14} M_{\odot} h^{-1}}, \qquad (\text{A.2})$$

while the specific values for the parameters γ and δ and their redshift dependence are shown in Table A.1.

Interestingly, the depletion of baryons in the Magneticum simulations is more pronounced than in TNG300, a difference that increases towards lower redshifts. At $z = 0$, a depletion of about 10 percent is observed for Magneticum even for the most massive clusters, while the baryon content in TNG300 clusters saturates to the cosmic value already for $M_{\text{vir}} \approx 2 \times 10^{14} M_{\odot} h^{-1}$.

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Baryon fraction inside the virial radius as a function of the virial mass. Top: Magneticum suite of simulations. Orange and blue symbols correspond to Box 2 and Box 1a, respectively. Bottom: TNG300 simulation. Different columns correspond to the redshifts $z \in \{0.0, 0.5, 1.0, 2.0\}$. The assumed cosmic baryon fraction for each simulation is depicted as the horizontal black line.</figcaption></figure>