EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 698 (June 2025) A&A, 698 (2025) A237 Full HTML
Open Access
Issue
A&A
Volume 698, June 2025
Article Number A237
Number of page(s) 17
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202452568
Published online 18 June 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

One of the key discoveries of the last century was that the Universe is expanding. At present, we have a well-constrained standard model, Lambda Cold Dark Matter (ΛCDM), which states that the Universe is composed of dark energy (ΩΛ∼0.69) and matter (Ωm=Ωb+Ωdm∼0.31), where Ωb represents the baryonic density and Ωdm denotes the dark matter density. From the density of matter, 15.7% corresponds to baryons, while ∼84.3% is attributed to dark matter (Planck Collaboration XIII 2016). This predominance of dark matter over visible matter significantly impacts the gravitational processes driving the formation of structures. Briefly, after the Big Bang, small quantum perturbations in the density field led to matter overdensities. As dark matter only interacts gravitationally, it had time to clump up and grow, forming the first dark matter haloes. These regions started to grow by accreting baryons and smaller haloes (White & Rees 1978; White & Frenk 1991). Several authors have thoroughly reviewed this hierarchical structure formation over the last decade, through analytical studies (e.g. Press & Schechter 1974; White & Rees 1978; White & Frenk 1991) and numerical simulations (e.g. White 1978; Villumsen 1982; McGee et al. 2009; Pallero et al. 2019).

thumbnail Fig. 1.

Distribution of the percentage difference between the stellar masses from Subfind and 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2}\sim r_{\mathrm {op}} $ apertures under different constrains. The black, red, and blue lines correspond to measurements within halocentric radius R500c, R200c, and the mass considered by the FoF algorithm, respectively. The dashed lines represent the normal distribution fitted to each constraint and are shown in their corresponding colour.

Within these haloes, galaxy clusters play a key role in helping us to understand the evolution of the Universe, as they represent the largest gravitationally bound structures. To form this kind of structure, matter in large regions of the Universe collapse (∼10 Mpc), resulting in a baryon fraction close to the universal value f bar U ∼ Ω b / Ω m ∼ 0.157 ± 0.004 $ f_{\mathrm {bar}}^{\mathrm {U}} \sim \Omega _{\mathrm {b}}/\Omega _{\mathrm {m}} \sim 0.157 \pm 0.004 $ (Planck Collaboration XIII 2016). Observational results indicate that massive clusters ( log 10 ( M 200 c / M ⊙ ) ≳ 14.5 $ \log _{10}{\left (M_{\mathrm {200c}}/ {\rm M}_{\odot }\right )} \gtrsim 14.5 $) typically exhibit a baryonic fraction close to this universal value, consistent with expectations (Giodini et al. 2009). Furthermore, Sunyaev–Zel’dovich (SZ) surveys and X-ray measurements confirm similar results (Chiu et al. 2018, hereafter C18). However, low-mass clusters and groups of galaxies deviate from these expectations (e.g. Ettori 2003; McGaugh et al. 2010; Lovisari et al. 2015; Chiu et al. 2018).

The last point highlights a discrepancy between early-time observations of baryon fractions in the cosmic microwave background (CMB) (e.g. Planck Collaboration XIII 2016) and late-time SZ and X-ray observations of galaxy groups and clusters. For example, C18 analysed the baryon fraction, fbar,500c=fstar,500c+fhot−gas,500c, of 91 self-similar clusters in the redshift range 0.20<z<1.25, using data from the SPT-SZ survey and Chandra X-ray observations. Using scaling relations for X-ray hot gas and galaxy stellar mass, they found that clusters with log 10 ( M 500 c / M ⊙ ) ≳ 14.9 $ \log _{10}{\left (M_{\mathrm {500c}}/ {\rm M}_{\odot }\right )} \gtrsim 14.9 $ reach a baryon fraction close to the universal value. However, clusters with log 10 ( M 500 c / M ⊙ ) ≲ 14.9 $ \log _{10}{\left (M_{\mathrm {500c}}/ {\rm M}_{\odot }\right )} \lesssim 14.9 $ exhibit significantly lower baryon fractions. These observations suggest that not all baryonic matter in low-mass clusters is detected, even though it should reside within their halo boundaries (Shull et al. 2012; Sanderson et al. 2013). This discrepancy presents a key challenge in the ΛCDM model, commonly termed the missing baryon problem. Additionally, C18 found that baryon fractions vary strongly with total mass but weakly with redshift. These two trends are inconsistent with a hierarchical scenario, as the most massive clusters with baryon fractions near the universal value could not have formed through the accretion of low-mass clusters with lower baryon fractions.

Studies from cosmological simulations and observations led to two main explanations for the missing baryons. The first suggests the existence of warm-hot gas in the intergalactic medium (WHIM) at temperatures of 105−107 K (Cen & Ostriker 1999). Nicastro et al. (2018) detected diffuse WHIM using two OVII absorption line tracers, although their measurements involve considerable uncertainties, as the signals could originate from sources other than the WHIM itself (Johnson et al. 2019). Furthermore, de Graaff et al. (2019) detected WHIM through the Planck thermal SZ (tSZ) signal (Planck Collaboration XXII 2016), and WHIM has been characterised and analysed in magnetohydrodynamical simulations (e.g. Martizzi et al. 2019). The second explanation posits that the missing baryons reside in regions outside the halo boundary, such as R500c or R200c (e.g. Haider et al. 2016; Ayromlou et al. 2021, 2023; Gouin et al. 2022). Feedback processes redistribute baryons to the outskirts of the halo, particularly in haloes with log 10 ( M 200 c / M ⊙ ) ≲ 14 $ \log _{10}\left (M_{\mathrm {200c}} / M_{\odot } \right ) \lesssim 14 $ (Ayromlou et al. 2023, hereafter A23). Using simulations, A23 introduced the closure radius, Rc, as the radius where the baryon fraction approaches the universal value. For example, for haloes with log 10 ( M 200 c / M ⊙ ) ≲ 14 $ \log _{10}\left (M_{\mathrm {200c}} / M_{\odot } \right ) \lesssim 14 $, active galactic nucleus (AGN) feedback processes redistribute baryons to Rc∼1.5−2.5R200c (Ayromlou et al. 2023). These two explanations are not mutually exclusive and provide complementary insights into galaxy evolution and large-scale structure cosmology (Walker et al. 2019; Nicastro et al. 2022).

Scaling relations are powerful tools for investigating the power laws associated with observables in haloes. The self-similar model (Kaiser 1986) provides a mathematical framework for describing scaling relations in hydrostatic equilibrium structures. These tools enable the use of large samples of galaxy clusters to constrain cosmological parameters (Haiman et al. 2001; Holder et al. 2001; Carlstrom et al. 2002). Scaling relations have been extensively studied through X-ray observations (e.g. Mohr & Evrard 1997; Arnaud & Evrard 1999; Mohr et al. 1999; Reiprich & Böhringer 2002; O’Hara et al. 2006; Arnaud et al. 2007; Pratt et al. 2009; Sun et al. 2009; Vikhlinin et al. 2009; Mantz et al. 2016; Chiu et al. 2018) and hydrodynamical simulations (e.g. Evrard 1997; Bryan & Norman 1998; Nagai et al. 2007; Stanek et al. 2010; Pillepich et al. 2018c; Pop et al. 2022a; Hadzhiyska et al. 2023a). The inclusion of AGN feedback processes in simulations significantly affects the X-ray and SZ scaling relations in galaxy groups (Puchwein et al. 2008; Fabjan et al. 2011; Pike et al. 2014; Planelles et al. 2013; Le Brun et al. 2016; Truong et al. 2017; Henden et al. 2018, 2019; Lim et al. 2021; Yang et al. 2022). Pop et al. (2022a, hereafter P22) introduced a novel method for fitting a scaling relation using a general double power-law (GDPL) model derived from the BPL (Jóhannesson et al. 2006) scheme. The GDPL model demonstrated superior performance in representing self-similarity for both relaxed and unrelaxed haloes at log 10 ( M 500 c / M ⊙ ) ≳ 14 $ \log _{10}\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right ) \gtrsim 14 $ compared to the simple power law (SPL: Lotka 1926) in simulations. If the trends and slopes of the SPL scaling relations presented by C18 can be reproduced in simulations, the missing baryons can be identified and further investigated.

In this study, we compared the self-similar scaling relations of galaxy clusters obtained from the observations of C18 with the hydrodynamical simulations. We used the suite of cosmological magnetohydrodynamical simulations from the IllustrisTNG project, selecting a sample of haloes with M200c≥7×1013 M⊙ to match the observational mass conditions and scaling relations of C18. To study their behaviour separately, we introduced cold and warm gas components into these scaling relations. We analysed the percentage distribution of the baryonic components of the haloes: cold gas, WHIM, hot gas, stellar mass from galaxies, and intra-cluster light (ICL). Finally, we fit the SPL, BPL, and GDPL models and examine their respective shapes.

This paper is structured as follows. In Section 2, we describe the methods used in this research, including simulations and observational data. Section 3 outlines the formalism of baryonic self-similar scaling relations, incorporating both the SPL and GDPL approaches. In Section 4, we present the results of the SPL, BPL, and GDPL models for the baryonic fraction scaling relations, comparing these models with observational data, varying the galaxy aperture, and adding the remaining gas components. In Section 5, we discuss the implications of these results, particularly in the context of baryonic component accretion. Finally, in Section 6, we summarise our findings.

2. Simulations and observational data

2.1. IllustrisTNG: The state-of-the-art magnetohydrodynamical and cosmological simulations

IllustrisTNG (Naiman et al. 2018; Nelson et al. 2018; Springel et al. 2018; Pillepich et al. 2018c; Marinacci et al. 2018) is a comprehensive suite of state-of-the-art magnetohydrodynamical and cosmological simulations focused on the formation and evolution of galaxies. It uses the moving-mesh AREPO code (Springel 2010), which solves the Poisson equations to compute the gravitational potential. In addition, it solves magnetohydrodynamic equations using a Voronoi tessellation method (Pakmor et al. 2011; Pakmor & Springel 2013). It includes a physical model for galaxy formation, incorporating non-gravitational processes such as gas radiative cooling, star formation, stellar evolution, and stellar feedback (Pillepich et al. 2018b). IllustrisTNG also accounts for the physics of supermassive black holes, including their seeding, merging, and feedback processes (Weinberger et al. 2017).

This work uses the TNG100 (∼100 Mpc box size) and TNG300 (∼300 Mpc box size) simulations. TNG100 offers a good compromise between a wide range of environments and resolution, while TNG300 provides the lowest resolution but allows for superior statistical analyses of galaxy clusters (see Table 1).

Table 1.

Parameters of TNG100 and TNG300 simulations from IllustrisTNG (Pillepich et al. 2018c).

We define galaxy clusters as self-gravitating haloes with masses M200c≥7×1013 M⊙ (for any redshift), comparable to the mass of the Fornax cluster (Drinkwater et al. 2001). The haloes and subhaloes were identified in a two-step fashion. First, haloes were selected using the FoF algorithm (FoF: Davis & Djorgovski 1985; Springel et al. 2001), with a linking length of b = 0.2, applied only to dark matter particles. Then, subhaloes were selected within FoF haloes using the Subfind algorithm (Springel et al. 2001; Dolag et al. 2009), as those self-gravitating structures residing within haloes. In this work, we define galaxies as those subhaloes containing 1000 stellar particles, corresponding to Mstar≥5×109 M⊙ for TNG100 and Mstar≥5×1010 M⊙ for TNG300, to ensure sufficient resolution to avoid spurious identifications. The gas within the structures was divided into three components: cold gas (T<105 K), warm gas (105≤T<107 K), and hot gas (T≥107 K) (Cen & Ostriker 1999; Martizzi et al. 2019; Gouin et al. 2022, 2023). These temperature thresholds were chosen because the star formation model in IllustrisTNG uses an effective equation of state and a representation of cold and warm gas phases within each cell (Pillepich et al. 2018b). Moreover, these thresholds are widely employed in IllustrisTNG studies to distinguish gas phases based on temperature and electron density (e.g. Martizzi et al. 2019; Gouin et al. 2022, 2023), as they effectively capture the gas phase diagram. Clusters were selected within 0.0≤z≤1.5, with a time window between snapshots of 0.5 Gyr (TNG100) and 1.8 Gyr (TNG300). This resulted in a sample of 218 and 1605 clusters, respectively.

For each cluster, we defined a radius RΔ=R500c, where the density is 500 times the critical density of the Universe (see Section 3). The baryon (gas and stars), stellar mass, and the hot, warm and cold gas masses were defined by summing the corresponding particles with RΔ. The total mass, M500c, is the sum of Mbar,500c and the dark matter mass, Mdm,500c. Similarly, we used the cluster total mass, which was determined by the FoF algorithm. These masses were defined as Mbar,FoF for baryons, Mdm,FoF for dark matter, and MFoF≡Mbar,FoF+Mdm,FoF for the total mass. The masses measured within R500c were used to compare the scaling relations in the simulations with the observations, while the FoF masses were used to examine whether the self-similarity of the scaling relations was satisfied.

For galaxies in clusters, we defined two apertures for mass measurement. The first corresponded to the Subfind-measured masses of stars, gas (hot, warm, cold, or total), and dark matter. The second corresponded to a spherical aperture with a radius R sph ≡ 2 × R M star / 2 $ R_{\mathrm {sph}} \equiv 2 \times R_{M_{\mathrm {star}}/2} $, comparable to the optical radius Rsph∼rop, which encloses most of the stellar flux in spectral energy distribution (SED) fitting (hereafter Rsph≡rop). This second aperture was beneficial for comparing simulations with observations and for estimating the impact of systematic uncertainties in the scaling relations due to differences in the metrics to measure masses. Finally, we characterised the intracluster medium (ICM) as all mass within the cluster not belonging to galaxies, as per our definitions. Table 2 summarises all mass notations used in the present study.

Table 2.

Mass notations.

2.2. Galaxy clusters in observations

To compare our results, we used the observations of C18, which analysed 91 galaxy clusters identified via the SZ effect (Sunyaev & Zeldovich 1972) in the redshift range 0.25≤z≤1.25. This sample was selected from the SPT-SZ survey (Bleem et al. 2015) and matched with Chandra X-ray observations. These measurements were further correlated with optical photometry in the griz bands from the Dark Energy Survey (DES; Dark Energy Survey Collaboration 2016) and near-infrared photometry obtained using either the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) or the Infrared Array Camera (IRAC; Fazio et al. 2004) on the Spitzer telescope.

C18 used the ζ−M500 scaling relations for the SZ effect observable, provided by the SPT collaboration, and the ζ normal distribution with the signal-to-noise ratio of the SZ effect ξ as the unit width. The calibrations of de Haan et al. (2016) were used to estimate the total mass, M500c, for each cluster in the sample (based on Bocquet et al. 2015). The mass range of these clusters was M 500 c = ( 2.63 − 12.71 ) × 10 14 M ⊙ $ M_{\mathrm {500c}} = \left (2.63 {-} 12.71 \right ) \times 10^{14}\,{\rm M}_{\odot } $. The hot-gas mass, Mhot−gas,500c, was determined using Chandra X-ray observations by fitting the best-fitting modified β-model (Vikhlinin et al. 2009) out to R500c (based on McDonald et al. 2013). The stellar mass, Mstar,500c, was estimated by fitting the SEDs of individual galaxies, constructing a stellar mass function (SMF; Schechter 1976), and summing the stellar masses of all galaxies in the cluster. Finally, the baryon fractions were calculated as the sum of the X-ray gas mass and the stellar mass of galaxies, normalised by M500c. The baryon fractions for these clusters ranged from fChiu,500c≈0.07−0.21.

3. Scaling relations

This section mainly uses the mathematical frameworks of Kaiser (1986), Mulroy et al. (2019), and P22. We followed their mathematical framework for the self-similar scaling relations of baryonic fractions in Section 3.1. Finally, the fitting procedures for the scaling relations are defined in Section 3.2.

3.1. Self-similar model and baryonic mass fraction scaling relations

We assumed that the hot gas was in hydrostatic equilibrium. Thus, we expected the pressure, Pgas, and the density, ρgas, of the gas to satisfy Eq. (1):

1 ρ gas d P gas d r = − GM r 2 , $$ \dfrac {1}{\rho _{\mathrm {gas}}} \dfrac {{\mathrm {d}} P_{\mathrm {gas}}}{{\mathrm {d}}r} = - \dfrac {GM}{r^{2}}, $$(1)

where the solution to the previous equation for the total mass, as a function of the gas density ρgas(r) and temperature T(r) profiles within the cluster, is given in Eq. (2):

M ( < r ) = − k B rT ( r ) μ m p G ( d ln ρ gas ( r ) d ln r + d ln T ( r ) d ln r ) . $$ M( < r) = -\dfrac {k_{\mathrm {B}} r T(r)}{\mu m_{\mathrm {p}} G} \left (\dfrac {{\mathrm {d \ln }} \rho _{\mathrm {gas}}(r)}{{\mathrm {d \ln }}r} + \dfrac {{\mathrm {d \ln }} T(r)}{{\mathrm {d \ln }} r} \right ). $$(2)

If the mass is defined on the basis of a spherical collapse model, the mass enclosed within a spherical region of radius RΔ, with Δ times the critical density of the Universe ρc(z) can be derived. This is given by Eq. (3):

M Δ = 4 π 3 Δ ρ c ( z ) R Δ 3 . $$ M_{\Delta } = \dfrac {4\pi }{3} \Delta \rho _{\mathrm {c}}(z) R_{\Delta }^{3}. $$(3)

Eq. (4) defines the critical density of the Universe, where H(z) = H0E(z), with H0 being the Hubble constant and E ( z ) ≡ Ω m ( 1 + z ) 3 + Ω Λ $ E(z) \equiv \sqrt {\Omega _{\mathrm {m}} \left (1 + z \right )^{3} + \Omega _{\Lambda }} $.

ρ c ( z ) = 3 H ( z ) 2 8 π G . $$ \rho _{\mathrm {c}}(z) = \dfrac {3 H(z)^{2}}{8\pi G}. $$(4)

It is important to mention that the definitions of RΔ and MΔ depend on the overdensity value Δ. In this study, we considered Δ={500,200}.

In the self-similar model (Kaiser 1986), the properties of baryonic matter are scaled solely by the total cluster mass, independent of redshift. Thus, Mstar,Δ∝Mgas,Δ∝Mstar,Δ+Mgas,Δ=Mbar,Δ∝MΔ. We defined the fraction of a given baryonic component X (stellar, gas, or baryonic = stellar + gas) as shown in Eq. (5):

f X , Δ ≡ M X , Δ / M Δ . $$ f_{X,\Delta } \equiv M_{X,\Delta } / M_{\Delta }. $$(5)

The scaling relation for these fractions is given in Eq. (6):

f star , Δ ∝ f gas , Δ ∝ f bar , Δ ∝ M Δ 0 = 1 , $$ f_{\mathrm {star,\Delta }} \propto f_{\mathrm {gas,\Delta }} \propto f_{\mathrm {bar,\Delta }} \propto M_{\Delta }^{0} = 1, $$(6)

however, it is important to understand that the self-similarity model can be disrupted by non-gravitational processes such as gas cooling, stripping, winds, and feedback (both stellar and AGN) (Kaiser 1986).

3.2. Fitting the power law to scaling relations

3.2.1. Simple power law

The self-similarity model assumes that the scaling relations of galaxy clusters follow a simple power-law (SPL) distribution (SPL: Lotka 1926) (Kaiser 1986). If we have two observables, (X,Y), the SPL model relates them using two free parameters, (ASPL,αSPL). Eq. (7) presents the simple power-law model:

Y = A SPL X α SPL . $$ Y = A_{\mathrm {SPL}} X^{\alpha _{\mathrm {SPL}}}. $$(7)

3.2.2. Broken power law and general double power law

One of the findings of P22 was that an SPL is insufficient to describe the scaling relations in IllustrisTNG simulations. Feedback processes are more intense in low-mass haloes and groups of galaxies (Ayromlou et al. 2021), breaking the self-similarity of the scaling relations for these low-mass haloes. To address this, we used the general double power-law (GDPL) model. This model is based on the broken power law (BPL: Jóhannesson et al. 2006), in which the data are divided into two regions, separated by a pivot point Xpivot. This value determines the slope values: αBPL,1 for X≤Xpivot and αBPL,2 for X>Xpivot. In the GDPL model, the transition between the slopes αGDPL,1 and αGDPL,2 is smoothed by a parameter δ, which quantifies the width of the transition region between the slopes (Pop et al. 2022a).

The BPL model was defined as

Y = A BPL · { ( X X pivot ) α BPL , 1 X ≤ X pivot ( X X pivot ) α BPL , 2 X > X pivot , $$ Y = A_{\mathrm {BPL}} \cdot \begin {cases}\left (\frac {X}{X_{\mathrm {pivot}}}\right )^{\alpha _{\mathrm {BPL,1}}} & X\leq X_{\mathrm {pivot}} \\ \left (\frac {X}{X_{\mathrm {pivot}}}\right )^{\alpha _{\mathrm {BPL,2}}} & X >X_{\mathrm {pivot}} \end {cases}, $$(8)

where the inclusion of a factor X pivot − 1 $ X_{\mathrm {pivot}}^{-1} $ in both pieces ensures that the two branches coincide at X=Xpivot, thereby guaranteeing continuity at the pivot point. The GDPL model was defined as

Y = A GDPL ( X X pivot ) α GDPL , 1 { 1 2 [ 1 + ( X X pivot ) 1 δ ] } ( α GDPL , 2 − α GDPL , 1 ) δ , $$ Y = A_{\mathrm {GDPL}} \left (\frac {X}{X_{\mathrm {pivot}}} \right )^{\alpha _{\mathrm {GDPL,1}}} \left\{ \frac {1}{2} \left [1 + \left (\frac {X}{X_{\mathrm {pivot}}} \right )^{\frac {1}{\delta }} \right ] \right\} ^{\left (\alpha _{\mathrm {GDPL,2}} - \alpha _{\mathrm {GDPL,1}}\right ) \delta }, $$(9)

with the main advantage of the GDPL model over the SPL model being that it allowed us to identify where self-similarity in the scaling relations holds. The parameter Xpivot defined the breakpoint between the first and second slopes, delineating the mass range where self-similarity breaks. Although introducing additional slopes with multiple breakpoints could improve the model's fitting performance, we found that a single Xpivot with two slopes was sufficient to capture the effects of AGN feedback within the mass range of our sample, as demonstrated in P22.

4. Results

4.1. Comparison between observations and IllustrisTNG

C18 measured halo masses using the SZ-effect, hot gas using the X-ray β-model, and the stellar masses of each galaxy through SED fitting and SMF construction, all within R500c. Thus, to enable a direct comparison, we considered two measurement approaches within the simulations. First, we included only the gas from the cluster with T≥107 K enclosed within a sphere of radius R500c. The stellar mass in this case was the sum of the stellar masses of galaxies within a spherical aperture of twice the stellar half-mass radius ( 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $). In this approach, warm gas, cold gas, and the ICL were not considered. Second, we considered all gas and stars within the halo inside a sphere of radius R500c. This included hot gas, warm gas, cold gas, galaxy stars, and the ICL.

The first approach was essential for comparing TNG100 and TNG300 with the observations of C18. The second approach enabled examination of the scaling relation behaviour when all components expected to be in the clusters are included. We did not use the Subfind aperture, as our goal was to create a measurement procedure that closely matches the observations and to characterise all baryonic components within the clusters. Fig. 1 shows the differences between the stellar masses defined by Subfind and 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ apertures within R500c, R200c, and FoF constraints. For TNG100, the difference was approximately 31.88±1.94 per cent, and for TNG300, 32.53±2.74 per cent within R500c, comparable to values observed for other constraints. These differences are significant for calculating baryon fractions, fitting the SPL scaling relation, and comparing with the scaling relations of C18. For masses within R200c and FoF, the percentages remained unchanged relative to R500c, with values ranging from 30.92% and 32.21% in both simulations.

Next, we tested the self-similarity of IllustrisTNG haloes. This was intended to remove the redshift dependence from the baryon fraction, allowing us to fit an SPL that depended solely on the mass, as described by Eq. (6). We calculated the baryon fraction using Eq. (5) for TNG100 and TNG300 haloes within a R500 as described above. The SPL used in this case was the relation f bar , 500 c ∝ ( 1 + z ) α SPL $ f_{\mathrm {bar, 500c}} \propto \left (1 + z\right )^{\alpha _{\mathrm {SPL}}} $. Fig. 2 shows the SPL fits of the IllustrisTNG haloes for the relation between baryon fraction and redshift, considering both comparison types. In the first comparison, the slopes were αSPL≃0.029±0.051 for TNG100 and αSPL≃0.071±0.021 for TNG300. In the second type of comparison, the slopes were αSPL≃0.006±0.023 for TNG100 and αSPL≃0.037±0.010 for TNG300. These slopes were obtained using a non-linear least squares regression fit to our SPL model. The errors on the slopes were determined from the covariance matrix of the fit. In the self-similar model presented by Kaiser (1986), the theoretical slope should be close to zero, as the baryonic matter should only scale with the total cluster mass, independently of redshift. Upon examining the slope values closely, we found that all were positive, indicating a slight change in the scaling relations, with a minimal increase in the baryon fraction with redshift. However, the differences between the self-similar and simulation slopes were small, on the order of hundredths (and in one case, thousandths). We therefore confirmed that the simulation slopes were effectively zero, removing the redshift dependence from the baryon fraction. As a result, the baryon fraction fbar,500c was expected to scale only with the mass M500c.

thumbnail Fig. 2.

SPL fit of baryon mass fraction fbar in IllustrisTNG at R500c and redshift z. Black dots show individual haloes; the green region indicates the SPL fit. Each panel is labelled with their corresponding simulation tag and the value for the αSPL slope from the SPL fit in the upper-left corner. Top panels: The fraction of baryons includes the hot gas and the galaxy stellar mass within the 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ aperture. Bottom panels: fraction of baryons includes all the gas and all stars inside the halo. The shaded areas are represented with a width of 3σ.

Using the same baryon fractions, we performed an SPL fit for IllustrisTNG haloes to determine the slopes, αSPL, and compared them with the data from C18. The SPL used in this case is the relation f bar , 500 c ∝ ( M 500 c ) α SPL $ f_{\mathrm {bar, 500c}} \propto \left (M_{\mathrm {500c}}\right )^{\alpha _{\mathrm {SPL}}} $. Fig. 3 presents these results along with the baryon fraction of the universe measured by Planck Collaboration XIII (2016). The top panels show that the observational data do not directly match TNG100 or TNG300, likely due to the AGN feedback model in IllustrisTNG, which efficiently enriches the intra-cluster medium gas of low-massive clusters and galaxy groups (Pop et al. 2022a). As a consequence, the baryon fraction increases at low mass values. However, this does not impact our analysis, as we focused on the behaviour of the scaling relation. We found that in IllustrisTNG, the relations between the baryon fraction and halo mass yielded a slope αSPL≃0.248±0.018 for TNG100 and αSPL≃0.261±0.006 for TNG300. Without considering redshift, C18 data provides an SPL scaling relation with a slope of αChiu≃0.350±0.076. The scaling relations for TNG100, TNG300, and C18 data show similar behaviour, with baryon fractions increasing at high halo masses and decreasing at low halo masses. Furthermore, the slope values of IllustrisTNG and C18 agree within 1−2σ. This agreement suggests that the IllustrisTNG results may provide insights into the missing baryons problem and help reconcile the tension with the hierarchical scenario highlighted by C18.

thumbnail Fig. 3.

SPL fit of baryon mass fraction fbar in IllustrisTNG at R500c for all redshifts. Black dots show individual haloes; the green region indicates the SPL fit. Each panel is labelled with their corresponding simulation tag and the value for the αSPL slope from the SPL fit in the upper-left corner. Top panels: fraction of baryons includes the hot gas and the galaxy stellar mass within the 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ aperture. Bottom panels: fraction of baryons includes all the gas and all the stars inside the halo. The yellow areas represent the Planck Collaboration XIII (2016) fraction of baryons. The C18 data (blue dots) is also shown, with its SPL fit (blue areas). The shaded areas are represented with a width of 3σ.

In the bottom panels of Fig. 3, we show how the behaviour of the scaling relations changes when the remaining gas (T<107 K) and all remaining stellar mass in the halo are included. A direct comparison with C18 is not possible, as the parameters are not directly comparable. Although the general behaviour remains, showing an increase in the baryon fraction, the slopes differ significantly between simulations and observations. We found a slope αSPL≃0.071±0.010 for TNG100 and αSPL≃0.120±0.003 for TNG300 in the scaling relations. These lower slopes result from including the remaining gas and accounting for the ICL in the measurements, which benefits low-mass haloes more than high-mass ones. Consequently, low-mass haloes exhibit a greater fraction of cold gas compared to high-mass haloes (see Section 4.2 and Appendix A). Additionally, more haloes in both simulations reach the universal baryon fraction. Most haloes that do not reach this value have masses of log10(M500c/M⊙)≲14. This mass threshold aligns with values reported in numerous studies discussing the breakdown of self-similarity in scaling relations (Pop et al. 2022a) and the redistribution of baryons in low-mass haloes by non-gravitational processes, such as AGN feedback (Ayromlou et al. 2023).

4.2. Looking for the missing baryons

A common limitation in observations is the lack of information on the cold and warm gas components. As simulations do not have this limitation, we estimated the fraction of baryons that remain undetected when only the hot gas in the ICM is considered. We divided the baryonic matter into cold, warm, and hot gas for the gaseous component, and into galaxy stellar mass (from rop and with Subfind) and the ICL for the stellar component. We further separated our clusters into low-mass log10(M500c/M⊙)<14 and high-mass haloes log10(M500c/M⊙)≥14. Our sample comprises 137 and 1115 low-mass haloes in IllustrisTNG-100 and IllustrisTNG-300, and 81 and 490 high-mass haloes in IllustrisTNG-100 and IllustrisTNG-300, respectively. For visualisation purposes, we also include the complete sample (labelled “all halo types”) in the following figures.

Table 3 presents the baryonic mass fraction for each component in our clusters. The baryonic mass fraction is defined as the mass of each component (Mτ), divided by the total baryonic mass (Mbar), where τ={cold gas, warm gas, hot gas, galaxy stars, ICL}. (see Appendix A and Figs. A.1. and A.2. for calculation and fitting details).

Table 3.

Percentages of baryonic components.

Furthermore, Table 3 shows the average baryonic mass values for components not considered in the C18 observations: cold gas, warm gas, and ICL (rop). We defined the sum of those components as the missing baryons in simulations. The data in Table 3 reveal a tendency for low-mass haloes to contain more missing (undetected) baryons than high-mass haloes. Specifically, low-mass haloes exhibit missing baryon fractions of approximately 17.83 per cent in TNG100 and 17.67 per cent in TNG300. In contrast, high-mass haloes have missing baryon fractions of about 6.95 per cent in TNG100 and 5.90 per cent in TNG300. In both cases, hot gas remains the most significant baryonic component within the haloes. However, warm gas and ICL (rop) correspond to a non-negligible fraction of baryons missed in low-mass haloes. In high-mass haloes, most of the missing baryons are attributed to the ICL (rop), which shows a similar proportional contribution to that observed in low-mass haloes. Finally, Table 3 also presents the stellar components measured using the Subfind aperture, highlighting the differences between this method and the rop aperture in terms of stellar mass measurements. The Subfind algorithm collects more stellar mass within galaxies compared to the rop spherical aperture, as it assigns all particles within the virial radius that are not associated with satellite subhaloes to the central galaxy (Springel et al. 2001; Dolag et al. 2009). As a result, diffuse stellar components, such as the stellar halo, are included as part of the central galaxy mass.

When considering all haloes, Table 3 shows that TNG100 and TNG300 have missing baryon fractions of 13.79 per cent and 14.08 per cent, respectively. Hot gas remains the most significant baryonic component. The overall fractions are similar to those of low-mass haloes but with lower contributions from warm gas and the ICL (rop) components, as high-mass haloes are included in this analysis. Warm gas is the most affected component when considering all haloes, showing a reduction of approximately three per cent compared to low-mass haloes in both simulations. The ICL (rop) component also exhibits a reduction of about 0.40 per cent relative to the low-mass haloes in both simulations. This analysis highlights the importance of dividing the halo sample, as the baryonic content can differ between high- and low-mass clusters.

Gouin et al. (2022) studied the distribution of gas components in TNG300 haloes at z = 0. They categorised these haloes by mass into clusters (M200c/M⊙/h>1×1014) and groups (5×1013<M200c/M⊙/h<1×1014) of galaxies. Additionally, the gas was divided into five independent phases: hot gas, warm dense gas (WCGM), warm diffuse gas (WHIM), cold dense gas (halo gas), and cold diffuse gas (IGM). The gas fraction of each phase was calculated as a function of R200. Table 4 presents the relative gas percentages for our haloes, calculated from the percentages in Table 3. Despite the differences between our sample and that of Gouin et al. (2022), the gas fractions are remarkably similar at R∼0.6R200∼R500. In galaxy groups, most of the contribution to these percentages, apart from hot gas, comes from warm diffuse gas, which corresponds to the circumgalactic medium (CGM) (Gouin et al. 2022). At this radius, WCGM gas accounts for a larger fraction than WHIM gas. In galaxy clusters, the majority of the gas fraction comes from hot gas. This analysis suggests that the gas phase associated with the missing baryons in our sample is likely dominated by WCGM gas.

Table 4.

Gas mass percentages.

These results suggest that the break in the power law seen in Fig. 3 is driven by the higher fractions of ICL, cold gas, and warm gas present in low-mass clusters compared to massive ones. The percentages of cold gas are lower in all cases compared to those of warm gas. Therefore, warm gas and the ICL in low-mass haloes are the primary contributors to the observed changes in baryonic scaling relations. This interpretation aligns with the WHIM proposal, which suggests that the missing baryons identified in C18 could reside in the warm diffuse phase. Furthermore, A23 propose that missing baryons are redistributed by non-gravitational processes (e.g. AGN, stellar, and UV feedback) into a region between 1.5 and 2.5R200c, where the baryon fraction reaches the universal value. Both scenarios are compatible and suggest a coherent explanation for the missing baryons problem. These findings may have significant implications for future observational research related to the baryon distribution, mass measurements, and the calibration of scaling relations.

4.3. FoF masses in baryonic scaling relations

We explored changes in baryonic scaling relations when considering the total mass of haloes within a sphere radius R500c, as discussed in Sections 4.1 and 4.2. The R500c was particularly useful for measuring gas masses in observations, especially those using X-ray and the SZ effect. In this section, we switch from using masses defined within R500c to those defined by the FoF algorithm. This approach accounts for all the mass bound to the cluster without a fixed aperture and allowed us to compare how scaling relations change when using these broader mass definitions.

Fig. 4 presents the SPL scaling relation models for baryonic mass fractions of IllustrisTNG haloes in the FoF-delimited mass context. Comparing these results with the all-gas + all-stars scaling relations at R500c shown in Fig. 3, we observed a mean slope difference of αSPL,FoF−αSPL,500c≃−0.06 in TNG100 and αSPL,FoF−αSPL,500c≃0.08 in TNG300. In the FoF context, low-mass haloes increase their baryon fractions, displaying a tendency towards αSPL,FoF→0 slopes in both simulations. This behaviour accounts for the differences observed with the R500c mass context. Furthermore, the self-similarity of the baryon fraction is satisfied, as the αSPL,FoF→0 results are consistent with Eq. (6). Furthermore, these slope tendencies indicate that high-mass haloes can form through the assembly of low-mass haloes, reflecting the hierarchical formation scenario. However, most IllustrisTNG haloes do not reach the baryon fraction measured by Planck Collaboration XIII (2016) within the 3σPlanck region (with σPlanck = 0.004 being the error of the baryon fraction in Planck Collaboration XIII 2016). Most haloes exhibit baryon fractions approximately 0.1 dex below the mean value, with only a few exceptions. This result suggests the re-emergence of the missing baryon problem in the FoF context, which will be further discussed in Section 5.

thumbnail Fig. 4.

SPL fit (green region) of the IllustrisTNG haloes (black dots), to the scaling relation of the mass fraction of baryons fbar and halo mass M measured by the FoF algorithm, for all redshifts. Each panel is labelled with the corresponding simulation tag and the value of the αSPL slope in the upper-left corner. The baryon fraction includes all gas and stellar mass within each halo. The yellow region represents the baryon fraction from Planck Collaboration XIII (2016). All regions are in 3σ.

4.4. Fitting power-law models as scaling relations

Figure 3, shows the SPL model fits to the scaling relations of baryon fractions in IllustrisTNG haloes. This provides a framework for comparison with the scaling relations from C18. However, the distribution of scattered data and the SPL model fit show significant deviations. In the hot gas + rop star comparison, the data exhibit a knee-like shape with a break (or pivot) around the low- and high-mass halo limit. This feature arises due to AGN feedback, which redistributes baryons to the outskirts of low-mass haloes, reducing their baryon fractions and breaking self-similarity in these haloes (Ayromlou et al. 2023; Pop et al. 2022a). As a result, the SPL model is insufficient for studying scaling relations across a wide range of haloes in the context of self-similarity (Pop et al. 2022a). To address this limitation, models capable of identifying breakpoints in the baryonic scaling relations are needed to determine where self-similarity is preserved. In this section, we employ the BPL and GDPL models for this purpose.

The method we followed was similar to that used by P22. We performed a non-linear least squares regression to fit the SPL, BPL, and GDPL models to the median baryon fraction values in the IllustrisTNG (TNG100 and TNG300) samples. To compute the median values, we used 20 bins, as in P22, but applied them to a higher mass range (M500c∈[1012−2×1015] M⊙). From these fits, we obtained the power-law parameters αSPL, αBPL,1, αBPL,2, αGDPL,1, and αGDPL,2, along with the δ and Mpivot parameters from the GDPL model (see details in Equations (7), (8), and (9)). The parameter Mpivot represents a breakpoint, indicating whether self-similarity was maintained or broken.

Figure 5 shows the SPL, BPL, and GDPL fits to the median values of the TNG100 and TNG300 halo samples, including the Mpivot value and the baryonic mass fraction of the Universe (Planck Collaboration XIII 2016). All numerical values of the power-law parameters – αSPL, αBPL,1, αBPL,2, αGDPL,1, and αGDPL,2 – as well as the GDPL smoothness parameter δ, the Mpivot value, and the respective model χ2-test results are presented in Table 5. The mean values and their associated errors were obtained from the fits to the median values of the sample. To compute the χ2-test, we used the median values of the TNG100 and TNG300 halo samples along with the models fitted with their mean parameter values. As shown in Fig. 5, the BPL and GDPL models exhibit a shape that closely matches the distribution of IllustrisTNG haloes in the hot-gas + rop case. In contrast, the SPL model matches the distribution for intermediate masses well but overestimates the baryon fraction at both low- and high-mass extremes. In the case of all gas + all stars, all three models demonstrate similar shapes compared to the halo median values, although the SPL model slightly overestimates the baryon fraction at the low- and high-mass extremes, but this effect is less pronounced than in the hot-gas + rop case. The χ2-test results for the SPL model support these claims, showing significant differences compared to the BPL and GDPL models in the hot-gas + rop case, but no notable differences in the all-gas + all-stars case. As detailed in Table 5, the χ2 values for the BPL and GDPL models were very close, with differences of χ GDPL 2 − χ BPL 2 ∈ [ 0 , 0.001 ] $ \chi ^{2}_{\mathrm {GDPL}} - \chi ^{2}_{\mathrm {BPL}} \in [0, 0.001] $, which is effectively negligible for both models.

thumbnail Fig. 5.

SPL fit (green region), BPL fit (magenta region) and GDPL fit (red region) of the IllustrisTNG haloes (median of the data in black solid line), to the scaling relation of the mass fraction of baryons fbar and halo mass M in a radius R500c. Each panel is labelled with the corresponding simulation tag in addition to a vertical dark green dashed line that represents the Mpivot of the broken fits. The power values of αSPL SPL fits, αBPL,1 and αBPL,2 BLP fits, and αGDPL,1 and αGDPL,2 GDPL fits, with the GDPL smoothness parameter δ and Mpivot value are in the Table 5. Top panels: fraction of baryons includes the hot-gas and the 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ aperture. Bottom panels: baryon fraction includes all gas and stellar mass within each halo. The yellow region represent the Planck Collaboration XIII (2016) fraction of baryons. All the regions are in 3σ.

Table 5.

Fit values of SPL, BPL and GDPL models.

In Table 5, we find that the SPL parameter fits for TNG100 in both cases are closer to the parameters of the scaling relations in Fig. 3. However, the SPL parameter fits for TNG300 in both cases are less consistent with these scaling relation parameters. This discrepancy is due to the different fitting method used in this analysis. Since both simulations use the same number of bins, TNG300 has a larger binning mass range than TNG100 due to its broader halo mass range. We also find that the mean power-law parameters (α) for the BPL and GDPL models are identical, but the BPL model shows smaller deviations than the GDPL model. As shown in Fig. 5, these standard deviations are negligible within the IllustrisTNG halo mass range but become noticeable outside of this range. This effect is more apparent in TNG100 than in TNG300, as the median mass values in TNG300 exhibit a flatter shape, resulting in similar fits for the BPL and GDPL models. The Mpivot values in the BPL and GDPL models are identical, with mean values in the range log 10 ( M pivot / M ⊙ ) = 13.80 − 14.07 $ \log _{10}\left (M_{\mathrm {pivot}} /{\mathrm {M_{\odot }}} \right ) = 13.80 {-} 14.07 $. This range includes the log 10 ( M 500 c / M ⊙ ) ∼ 14 $ \log _{10}\left (M_{\mathrm {500c}} /{\mathrm {M_{\odot }}} \right ) \sim 14 $ mass limit reported by P22 as the threshold for the loss of self-similarity in scaling relations, and by A23 as the upper limit for baryon redistribution driven by AGN feedback. This strongly suggests that AGN feedback lowers the baryonic mass fraction in low-mass haloes. Despite considering all baryonic components within R500c, low-mass haloes still show reduced baryon fractions, while high-mass haloes approach the baryonic fraction value reported by Planck Collaboration XIII (2016) within 3σPlanck. This supports the hypothesis that warm gas, cold gas, and ICL are under-represented in observations. Finally, the mean softened transition parameter, δ, for the GDPL models was consistent across both simulations and comparison types, with δ∼0.25. The σδ values are not significant, as δ primarily governs the transition between α1 and α2.

The BPL and GDPL scaling relations provide improved characterisation of non-gravitational processes using the baryon fraction of haloes, particularly in the context of galaxy clusters. The SPL model slopes derived from IllustrisTNG are consistent with those reported by C18, supporting the analysis of similar observations and the comparison of slope results with simulations. In the observational domain, ongoing and future surveys such as the Dark Energy Survey (Dark Energy Survey Collaboration 2016), Euclid (Laureijs et al. 2011), the Vera C. Rubin Observatory's LSST (LSST Dark Energy Science Collaboration 2012), eROSITA (Merloni et al. 2012; Pillepich et al. 2012, 2018a; Predehl et al. 2021; Bulbul et al. 2022), Athena X-ray observatory (Nandra et al. 2013), SPT-3G (Benson et al. 2014; Bender et al. 2018), CMB-S4 (Abazajian et al. 2016), CMB-HD (Sehgal et al. 2019), and Advanced ACTpol (Henderson et al. 2016) could validate and benefit from these models by comparing them with current and future cosmological hydrodynamical simulations. For instance, upcoming projects such as TNG-Cluster simulations (Nelson et al. 2024; Ayromlou et al. 2024; Lee et al. 2024; Lehle et al. 2024; Rohr et al. 2024; Truong et al. 2024), the MillenniumTNG Project (Hernández-Aguayo et al. 2023; Pakmor et al. 2023; Barrera et al. 2023; Kannan et al. 2023; Hadzhiyska et al. 2023c, a, b; Bose et al. 2023; Contreras et al. 2023; Delgado et al. 2023; Ferlito et al. 2023), and the FLAMINGO simulations (Schaye et al. 2023; Kugel et al. 2023) could provide valuable insights. These simulations will also allow for the extension of this framework by incorporating a larger sample of haloes, enabling a deeper understanding of baryonic scaling relations and their implications for cosmology.

5. Discussions

5.1. Mock observations and the fraction of baryons

In observational studies, the measurement of masses is limited by the 2D projection of the sky, which introduces biases into the estimates. Among the methods used to estimate mass, the SZ effect is useful to calculate M500c, while the fitting of the X-ray model β is useful for calculating Mhot−gas,500c. In contrast, simulations often overestimate observable values when using 3D projections instead of 2D projections. However, simulations can generate mock observations that provide an estimate of the mean bias, b, between the mass computed in simulations, M 500 c SIM $ M_{\mathrm {500c}}^{\mathrm {SIM}} $, and the observable mass in mocks, M 500 c EST $ M_{\mathrm {500c}}^{\mathrm {EST}} $ (Pop et al. 2022b). The relationship between these quantities is defined as

b ≡ 1 − M 500 c EST M 500 c SIM . $$ b \equiv 1 - \frac {M_{\mathrm {500c}}^{\mathrm {EST}}}{M_{\mathrm {500c}}^{\mathrm {SIM}}}. $$(10)

Pop et al. (2022b) studied the biases associated with the SZ effect, bSZ, and X-ray techniques, bX, in mass measurements using the MOCK-X pipeline in IllustrisTNG (Barnes et al. 2021). Pop et al. (2022b) identified 30 000 simulated galaxy groups and clusters with M500 ranging from 1012 to 2×1015 M⊙ at z = 0. For haloes with M500>1013 M⊙, the median hydrostatic mass bias was found to be bSZ = 0.125±0.003, while the X-ray mass bias was bX = 0.170±0.004.

If we consider both values for the estimation of the total mass M500c, defined as M 500 c EST = ( 1 − b i ) M 500 c SIM $ M_{\mathrm {500c}}^{\mathrm {EST}} = (1 - b_{\mathrm {i}}) M_{\mathrm {500c}}^{\mathrm {SIM}} $, with bi={bX,bSZ}, we define f bar , 500 c EST ≡ M bar , 500 c / M 500 c EST $ f_{\mathrm {bar,500c}}^{\mathrm {EST}} \equiv M_{\mathrm {bar, 500c}} / M_{\mathrm {500c}}^{\mathrm {EST}} $ as the synthetic baryon fraction and f bar , 500 c SIM $ f_{\mathrm {bar,500c}}^{\mathrm {SIM}} $ as the simulated value. The relation between both baryon fractions can be written as f bar , 500 c EST ≡ ( 1 − b f ) f bar , 500 c SIM $ f_{\mathrm {bar,500c}}^{\mathrm {EST}} \equiv (1 - b_{\mathrm {f}})f_{\mathrm {bar,500c}}^{\mathrm {SIM}} $, where bf is the baryon fraction bias. The parameter bf can be related to bi, and its associated uncertainty σ b f $ \sigma _{b_{\mathrm {f}}} $, through the following expressions:

b f ≡ 1 − 1 1 − b i , $$ b_{\mathrm {f}} \equiv 1 - \frac {1}{1 - b_{\mathrm {i}}}, $$(11)

σ b f = σ b i ( 1 − b i ) 2 . $$ \sigma _{b_{\mathrm {f}}} = \frac {\sigma _{b_{\mathrm {i}}}}{\left (1 - b_{\mathrm {i}} \right )^2}. $$(12)

Using the values from Pop et al. (2022b), we find bf≈−0.1429±0.0039 (for bSZ) and bf≈−0.2048±0.0058 (for bX). Under these assumptions, we expect an overestimation in the baryon fraction from 14.29 to 20.48 per cent, which represents a significant departure from the results of C18. However, this overestimation could be reduced if the baryonic mass bias, bbar, could be estimated. A condition of bbar>bi would be required to produce a decrease in f bar , 500 c EST $ f_{\mathrm {bar,500c}}^{\mathrm {EST}} $ relative to f bar , 500 c SIM $ f_{\mathrm {bar,500c}}^{\mathrm {SIM}} $. Moreover, bbar could account for individual biases associated with the stellar mass and the various gas phases.

5.2. Mass accretion in FoF masses

The baryon fraction values derived from FoF in TNG100 and TNG300 differ from the Planck Collaboration XIII (2016) value by approximately ∼0.1 dex, as shown in Fig. 4. This corresponds to fTNG/(Ωb/Ωm)≃0.79, indicating that about ∼21% of baryons are missing within clusters compared to the universal value. The redistribution of baryons to the outskirts of haloes by AGN feedback could explain this discrepancy for low-mass clusters (Ayromlou et al. 2023). However, the issue persists at higher masses, where AGN feedback is insufficient to redistribute baryons to large distances. In this context, one might consider whether the relaxation state of clusters (relaxed or unrelaxed) contributes to this discrepancy. However, this seems unlikely, as the majority of haloes fall below the universal baryon fraction, regardless of their relaxation state.

The FoF clustering algorithm assigns all gravitationally bound particles to the same halo. Consequently, this approach does not recover the universal baryon fraction,. In contrast, methods that use large concentric radii can achieve this at 1.5−2.5R200c (Ayromlou et al. 2023). However, this method also presents challenges: it may include elements not gravitationally bound to the cluster and does not account for the cluster's shape, as it only considers elements within a sphere. It would be particularly interesting to investigate the relationship between the shape of galaxy clusters and their baryon fraction using the concentric radius method.

The differences between the baryon fraction reported by Planck Collaboration XIII (2016), and that found in IllustrisTNG haloes when using the FoF algorithm may be attributed to two main factors. First, galaxy clusters may be accreting dark matter efficiently. Second, they may be accreting gas inefficiently.

Both hypotheses are related to the mass accretion rate (MAR) of dark matter (MARdm) being higher than that of gas (MARgas). Pizzardo et al. (2023b) developed a caustic technique to measure the mass accretion rate of galaxy clusters using spherical shells. The radii of these spherical shells are determined by the radial velocity profile, where the radial velocity satisfies vrad≤A×vmin, with vmin being the minimum radial velocity in the cluster. The parameter A is an empirically determined adjustment parameter, chosen to accurately capture the infall region described causticly. Pizzardo et al. (2023b) selected A = 0.72 to have a r3c=MAR3D/MARc = 1.00±0.04, where MAR3D and MARC are mass accretion rates obtained by 3D projections and caustic description, respectively. Using TNG300, Pizzardo et al. (2023b) analysed 1697 FoF haloes with M200c≥1014 M⊙, covering a redshift range of 0.01≤z≤1.04, ultimately working with 1318 haloes. The remaining ∼22% of the haloes were excluded due to limitations in the reliable application of the caustic technique (Pizzardo et al. 2023a). A notable finding concerns the mass accretion rates of galaxies and dark matter within these haloes. The study found that the MARdms are approximately 6.5% lower than the galaxy mass accretion rates (MARgalaxies) (Pizzardo et al. 2023b). This discrepancy arises because the clustering amplitude of galaxies relative to dark matter is larger on smaller scales and at higher redshifts (Davis & Peebles 1983; Davis et al. 1985; Jenkins et al. 1998).

However, applying this model to gas presents additional challenges. In IllustrisTNG, the evolution of gas is tracked using a cell-based method (Pillepich et al. 2018c). With this approach, tracing the gas becomes challenging as its vectorial information is not retained. Consequently, a gas particle simulation suite is required for more accurate tracking. The FLAMINGO simulations (Schaye et al. 2023; Kugel et al. 2023), for example, use a gas particle method in simulation boxes of 1.0 and 2.8 Gpc. A comparison of these results with FLAMINGO simulations, including an analysis of gas dynamics and mass accretion rate will be addressed in a future study.

6. Summary and conclusions

In this study, we examined the scaling relation differences between the state-of-the-art magnetohydrodynamic and cosmological simulations from IllustrisTNG and the 91 galaxy clusters identified via the SZ effect in the SPT-SZ survey (Bleem et al. 2015), as reported by C18. Our simulation sample consisted of 218 haloes from TNG100 and 1604 haloes from TNG300, each with M200c≥7×1013 M⊙. The baryonic mass fraction was defined as a function of halo mass within a concentric sphere of radius R500, with galaxy components selected using the Subfind algorithm.

We defined two observational comparison cases: (1) incorporating only the hot gas within the halo and galaxy stars within a spherical aperture of 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $; and (2) including all baryonic components (gas and stars) in the halo, both analysed at R500c. To test the self-similarity of the haloes, we applied an SPL model and found that the slopes in both cases and simulations converge towards αSPL→0, indicating the elimination of redshift dependence in the scaling relations.

Our main findings can be summarised as follows.

  • The slopes of the SPL scaling relations (αSPL) in TNG100 and TNG300 haloes align with those reported by C18 within an error of 1−2σ for the “hot gas + rop stars” case. This suggests a potential solution to the missing baryons problem identified in C18 and the tensions within the hierarchical formation scenario.

  • When all baryonic components (hot gas, warm gas, cold gas, ICL, and galaxy stars) are included, the SPL scaling relations exhibit a null slope. Low-mass haloes ( log 10 ( M 500 c / M ⊙ ) ≲ 14 $ \log _{10} \left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right ) \lesssim 14 $) show significantly elevated baryonic fractions, but remain below the universal baryon fraction, whereas high-mass haloes ( log 10 ( M 500 c / M ⊙ ) ≳ 14 $ \log _{10} \left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right ) \gtrsim 14 $) achieve this value. This mass threshold correlates with the loss of self-similarity in the haloes and suggests that non-gravitational processes play a critical role in baryon redistribution.

  • The mean percentage of missing baryons is higher in low-mass haloes (∼17.67% in TNG100, ∼17.80 per cent in TNG300) compared to high-mass haloes (∼6.95 per cent in TNG100, ∼5.90 per cent in TNG300). This disparity accounts for the observed behavioural differences in scaling relations across different baryonic component configurations. Furthermore, the missing baryons in the IllustrisTNG haloes are primarily due to the exclusion of ICL and warm gas, particularly in low-mass haloes, whereas cold gas contributes negligibly. This finding strongly supports the warm-hot intergalactic medium (WHIM) scenario.

  • We find that discrepancies in galaxy apertures, particularly in the characterisation of ICL, can lead to variation of approximately five per cent in the inferred missing baryon fraction, depending on the stellar mass measurement model. This variability must be considered in observational studies.

We investigated the accretion of baryonic components within FoF haloes and found that approximately ∼21% of baryons are missing compared to the expected universal fraction (Planck Collaboration XIII 2016). While the redistribution of baryons to the outskirts of low-mass haloes (Gouin et al. 2022; Ayromlou et al. 2023) may account for this deficit, it does not explain the shortfall in high-mass haloes. We hypothesise that FoF haloes may either be efficiently accreting dark matter or inefficiently accreting gas. Although the dark matter mass accretion rate (MARdm) can be traced in IllustrisTNG using a caustic technique (Pizzardo et al. 2023a, b), the gas mass accretion rate (MARgas) cannot be similarly traced due to the absence of positional and vectorial information for gas cells. Future work will focus on following gas particles in SPH simulation suites, such as FLAMINGO, to track gas accretion more effectively and calculate MARgas.

Finally, we emphasised the importance of incorporating non-gravitational processes in scaling relation analyses, particularly through the BPL and GDPL models. These models enable the identification of critical mass breakpoints where self-similarity is either lost or preserved. Future observational surveys could validate and refine these models, particularly when compared with current and upcoming simulations of large haloes. Such efforts may expand this framework by integrating additional data, further improving our understanding of baryonic scaling relations and their implications for cosmology.

Data availability

The data from the IllustrisTNG simulations used in this work is publicly available on the IllustrisTNG website (https://www.tng-project.org/). The haloes were selected using the JupyterLab Workspace and their own IllustrisTNG modules. The data of the SPT-SZ survey clusters is available in Tables 1 and 2 of C18 work.

Acknowledgments

We thank the EVOLGAL4D Computational Galaxy Formation and Evolution team of Pontificia Universidad Católica de Chile for the support during all the research time (visit EVOLGAL4D website in https://www.evolgal4d.com/). We thank Facundo Gómez and Joseph J. Mohr for the early discussions about this topic. We thank the anonymous referee for their constructive comments and suggestions, which helped improve the clarity and quality of this manuscript. This work was supported by partial funding from the ANID Basal Project FB210003. DM acknowledges financial support from PUCV Maintenance and Tuition Grants. DP acknowledges financial support from ANID through FONDECYT Postdoctrorado Project 3230379. DP also gratefully acknowledges financial support from ANID – MILENIO – NCN2024_112. PBT acknowledges the partial funding by FONDECYT 1240465. This research made use of several software packages, including NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Pandas (McKinney et al. 2010), and all figures were generated using matplotlib (Hunter 2007).

References

  1. Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv e-prints [arXiv:1610.02743] [Google Scholar]
  2. Arnaud, M., & Evrard, A. E. 1999, MNRAS, 305, 631 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2007, A&A, 474, L37 [CrossRef] [EDP Sciences] [Google Scholar]
  4. Ayromlou, M., Nelson, D., Yates, R. M., et al. 2021, MNRAS, 502, 1051 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ayromlou, M., Nelson, D., & Pillepich, A. 2023, MNRAS, 524, 5391 [CrossRef] [Google Scholar]
  6. Ayromlou, M., Nelson, D., Pillepich, A., et al. 2024, A&A, 690, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Barnes, D. J., Vogelsberger, M., Pearce, F. A., et al. 2021, MNRAS, 506, 2533 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barrera, M., Springel, V., White, S. D. M., et al. 2023, MNRAS, 525, 6312 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bender, A. N., Ade, P. A. R., Ahmed, Z., et al. 2018, in Proc. SPIE, eds. J. Zmuidzinas, & J. R. Gao, 10708, 1070803 [Google Scholar]
  10. Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in Proc. SPIE, eds. W. S. Holland, & J. Zmuidzinas, 9153, 91531P [NASA ADS] [CrossRef] [Google Scholar]
  11. Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27 [Google Scholar]
  12. Bocquet, S., Saro, A., Mohr, J. J., et al. 2015, ApJ, 799, 214 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bose, S., Hadzhiyska, B., Barrera, M., et al. 2023, MNRAS, 524, 2579 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bulbul, E., Liu, A., Pasini, T., et al. 2022, A&A, 661, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Carlstrom, J. E., Holder, G. P., & Reese, E. D. 2002, ARA&A, 40, 643 [Google Scholar]
  17. Cen, R., & Ostriker, J. P. 1999, ApJ, 514, 1 [NASA ADS] [CrossRef] [Google Scholar]
  18. Chiu, I., Mohr, J. J., McDonald, M., et al. 2018, MNRAS, 478, 3072 [Google Scholar]
  19. Contreras, S., Angulo, R. E., Springel, V., et al. 2023, MNRAS, 524, 2489 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dark Energy Survey Collaboration (Abbott, T., et al.) 2016, MNRAS, 460, 1270 [Google Scholar]
  21. Davis, M., & Djorgovski, S. 1985, ApJ, 299, 15 [Google Scholar]
  22. Davis, M., & Peebles, P. J. E. 1983, ApJ, 267, 465 [Google Scholar]
  23. Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]
  24. de Graaff, A., Cai, Y. -C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. de Haan, T., Benson, B. A., Bleem, L. E., et al. 2016, ApJ, 832, 95 [NASA ADS] [CrossRef] [Google Scholar]
  26. Delgado, A. M., Hadzhiyska, B., Bose, S., et al. 2023, MNRAS, 523, 5899 [CrossRef] [Google Scholar]
  27. Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497 [Google Scholar]
  28. Drinkwater, M. J., Gregg, M. D., & Colless, M. 2001, ApJ, 548, L139 [Google Scholar]
  29. Ettori, S. 2003, MNRAS, 344, L13 [CrossRef] [Google Scholar]
  30. Evrard, A. E. 1997, MNRAS, 292, 289 [CrossRef] [Google Scholar]
  31. Fabjan, D., Borgani, S., Rasia, E., et al. 2011, MNRAS, 416, 801 [NASA ADS] [CrossRef] [Google Scholar]
  32. Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10 [Google Scholar]
  33. Ferlito, F., Springel, V., Davies, C. T., et al. 2023, MNRAS, 524, 5591 [CrossRef] [Google Scholar]
  34. Giodini, S., Pierini, D., Finoguenov, A., et al. 2009, ApJ, 703, 982 [NASA ADS] [CrossRef] [Google Scholar]
  35. Gouin, C., Gallo, S., & Aghanim, N. 2022, A&A, 664, A198 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Gouin, C., Bonamente, M., Galárraga-Espinosa, D., Walker, S., & Mirakhor, M. 2023, A&A, 680, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Hadzhiyska, B., Ferraro, S., Pakmor, R., et al. 2023a, MNRAS, 526, 369 [NASA ADS] [CrossRef] [Google Scholar]
  38. Hadzhiyska, B., Hernquist, L., Eisenstein, D., et al. 2023b, MNRAS, 524, 2524 [NASA ADS] [CrossRef] [Google Scholar]
  39. Hadzhiyska, B., Eisenstein, D., Hernquist, L., et al. 2023c, MNRAS, 524, 2507 [NASA ADS] [CrossRef] [Google Scholar]
  40. Haider, M., Steinhauser, D., Vogelsberger, M., et al. 2016, MNRAS, 457, 3024 [Google Scholar]
  41. Haiman, Z., Mohr, J. J., & Holder, G. P. 2001, ApJ, 553, 545 [NASA ADS] [CrossRef] [Google Scholar]
  42. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  43. Henden, N. A., Puchwein, E., Shen, S., & Sijacki, D. 2018, MNRAS, 479, 5385 [CrossRef] [Google Scholar]
  44. Henden, N. A., Puchwein, E., & Sijacki, D. 2019, MNRAS, 489, 2439 [NASA ADS] [CrossRef] [Google Scholar]
  45. Henderson, S. W., Allison, R., Austermann, J., et al. 2016, J. Low Temp. Phys., 184, 772 [NASA ADS] [CrossRef] [Google Scholar]
  46. Hernández-Aguayo, C., Springel, V., Pakmor, R., et al. 2023, MNRAS, 524, 2556 [CrossRef] [Google Scholar]
  47. Holder, G., Haiman, Z., & Mohr, J. J. 2001, ApJ, 560, L111 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hunter, J. D. 2007, Comp. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  49. Jenkins, A., Frenk, C. S., Pearce, F. R., et al. 1998, ApJ, 499, 20 [Google Scholar]
  50. Jóhannesson, G., Björnsson, G., & Gudmundsson, E. H. 2006, ApJ, 640, L5 [Google Scholar]
  51. Johnson, S. D., Mulchaey, J. S., Chen, H. -W., et al. 2019, ApJ, 884, L31 [NASA ADS] [CrossRef] [Google Scholar]
  52. Kaiser, N. 1986, MNRAS, 222, 323 [Google Scholar]
  53. Kannan, R., Springel, V., Hernquist, L., et al. 2023, MNRAS, 524, 2594 [CrossRef] [Google Scholar]
  54. Kugel, R., Schaye, J., Schaller, M., et al. 2023, MNRAS, 526, 6103 [NASA ADS] [CrossRef] [Google Scholar]
  55. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193] [Google Scholar]
  56. Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2016, MNRAS, 466, 4442 [Google Scholar]
  57. Lee, W., Pillepich, A., ZuHone, J., et al. 2024, A&A, 686, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Lehle, K., Nelson, D., Pillepich, A., Truong, N., & Rohr, E. 2024, A&A, 687, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Lim, S. H., Barnes, D., Vogelsberger, M., et al. 2021, MNRAS, 504, 5131 [CrossRef] [Google Scholar]
  60. Lotka, A. J. 1926, J. Wash. Acad. Sci., 16, 317 [Google Scholar]
  61. Lovisari, L., Reiprich, T. H., & Schellenberger, G. 2015, A&A, 573, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. LSST Dark Energy Science Collaboration 2012, ArXiv e-prints [arXiv:1211.0310] [Google Scholar]
  63. Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2016, MNRAS, 463, 3582 [NASA ADS] [CrossRef] [Google Scholar]
  64. Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
  65. Martizzi, D., Vogelsberger, M., Artale, M. C., et al. 2019, MNRAS, 486, 3766 [Google Scholar]
  66. McDonald, M., Benson, B. A., Vikhlinin, A., et al. 2013, ApJ, 774, 23 [NASA ADS] [CrossRef] [Google Scholar]
  67. McGaugh, S. S., Schombert, J. M., de Blok, W. J. G., & Zagursky, M. J. 2010, ApJ, 708, L14 [NASA ADS] [CrossRef] [Google Scholar]
  68. McGee, S. L., Balogh, M. L., Bower, R. G., Font, A. S., & McCarthy, I. G. 2009, MNRAS, 400, 937 [Google Scholar]
  69. McKinney, W. 2010, in Proc. 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56–61 [CrossRef] [Google Scholar]
  70. Merloni, A., Predehl, P., Becker, W., et al. 2012, ArXiv e-prints [arXiv:1209.3114] [Google Scholar]
  71. Mohr, J. J., & Evrard, A. E. 1997, ApJ, 491, 38 [Google Scholar]
  72. Mohr, J. J., Mathiesen, B., & Evrard, A. E. 1999, ApJ, 517, 627 [Google Scholar]
  73. Mulroy, S. L., Farahi, A., Evrard, A. E., et al. 2019, MNRAS, 484, 60 [NASA ADS] [CrossRef] [Google Scholar]
  74. Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1 [Google Scholar]
  75. Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
  76. Nandra, K., Barret, D., Barcons, X., et al. 2013, ArXiv e-prints [arXiv:1306.2307] [Google Scholar]
  77. Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 [Google Scholar]
  78. Nelson, D., Pillepich, A., Ayromlou, M., et al. 2024, A&A, 686, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Nicastro, F., Kaastra, J., Krongold, Y., et al. 2018, Nature, 558, 406 [Google Scholar]
  80. Nicastro, F., Fang, T., & Mathur, S. 2022, ArXiv e-prints [arXiv:2203.15666] [Google Scholar]
  81. O’Hara, T. B., Mohr, J. J., Bialek, J. J., & Evrard, A. E. 2006, ApJ, 639, 64 [Google Scholar]
  82. Pakmor, R., & Springel, V. 2013, MNRAS, 432, 176 [NASA ADS] [CrossRef] [Google Scholar]
  83. Pakmor, R., Bauer, A., & Springel, V. 2011, MNRAS, 418, 1392 [NASA ADS] [CrossRef] [Google Scholar]
  84. Pakmor, R., Springel, V., Coles, J. P., et al. 2023, MNRAS, 524, 2539 [NASA ADS] [CrossRef] [Google Scholar]
  85. Pallero, D., Gómez, F. A., Padilla, N. D., et al. 2019, MNRAS, 488, 847 [NASA ADS] [CrossRef] [Google Scholar]
  86. Pike, S. R., Kay, S. T., Newton, R. D. A., Thomas, P. A., & Jenkins, A. 2014, MNRAS, 445, 1774 [NASA ADS] [CrossRef] [Google Scholar]
  87. Pillepich, A., Porciani, C., & Reiprich, T. H. 2012, MNRAS, 422, 44 [Google Scholar]
  88. Pillepich, A., Reiprich, T. H., Porciani, C., Borm, K., & Merloni, A. 2018a, MNRAS, 481, 613 [Google Scholar]
  89. Pillepich, A., Springel, V., Nelson, D., et al. 2018b, MNRAS, 473, 4077 [Google Scholar]
  90. Pillepich, A., Nelson, D., Hernquist, L., et al. 2018c, MNRAS, 475, 648 [Google Scholar]
  91. Pizzardo, M., Geller, M. J., Kenyon, S. J., Damjanov, I., & Diaferio, A. 2023a, A&A, 675, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Pizzardo, M., Geller, M. J., Kenyon, S. J., Damjanov, I., & Diaferio, A. 2023b, A&A, 680, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Planelles, S., Borgani, S., Fabjan, D., et al. 2013, MNRAS, 438, 195 [Google Scholar]
  96. Pop, A. R., Hernquist, L., Nagai, D., et al. 2022a, MNRAS, submitted [arXiv:2205.11528] [Google Scholar]
  97. Pop, A. R., Hernquist, L., Nagai, D., et al. 2022b, MNRAS, submitted [arXiv:2205.11537] [Google Scholar]
  98. Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1 [EDP Sciences] [Google Scholar]
  100. Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [Google Scholar]
  101. Puchwein, E., Sijacki, D., & Springel, V. 2008, ApJ, 687, L53 [NASA ADS] [CrossRef] [Google Scholar]
  102. Reiprich, T. H., & Böhringer, H. 2002, ApJ, 567, 716 [Google Scholar]
  103. Rohr, E., Pillepich, A., Nelson, D., Ayromlou, M., & Zinger, E. 2024, A&A, 686, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Sanderson, A. J. R., O’Sullivan, E., Ponman, T. J., et al. 2013, MNRAS, 429, 3288 [NASA ADS] [CrossRef] [Google Scholar]
  105. Schaye, J., Kugel, R., Schaller, M., et al. 2023, MNRAS, 526, 4978 [NASA ADS] [CrossRef] [Google Scholar]
  106. Schechter, P. 1976, ApJ, 203, 297 [Google Scholar]
  107. Sehgal, N., Aiola, S., Akrami, Y., et al. 2019, BAAS, 51, 6 [NASA ADS] [Google Scholar]
  108. Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23 [NASA ADS] [CrossRef] [Google Scholar]
  109. Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]
  110. Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726 [Google Scholar]
  111. Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676 [Google Scholar]
  112. Stanek, R., Rasia, E., Evrard, A. E., Pearce, F., & Gazzola, L. 2010, ApJ, 715, 1508 [NASA ADS] [CrossRef] [Google Scholar]
  113. Sun, M., Voit, G. M., Donahue, M., et al. 2009, ApJ, 693, 1142 [NASA ADS] [CrossRef] [Google Scholar]
  114. Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173 [NASA ADS] [EDP Sciences] [Google Scholar]
  115. Truong, N., Rasia, E., Mazzotta, P., et al. 2017, MNRAS, 474, 4089 [Google Scholar]
  116. Truong, N., Pillepich, A., Nelson, D., et al. 2024, A&A, 686, A200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  117. Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009, ApJ, 692, 1033 [Google Scholar]
  118. Villumsen, J. V. 1982, MNRAS, 199, 493 [Google Scholar]
  119. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  120. Walker, S., Simionescu, A., Nagai, D., et al. 2019, Space Sci. Rev., 215, 7 [Google Scholar]
  121. Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291 [Google Scholar]
  122. White, S. D. M. 1978, MNRAS, 184, 185 [NASA ADS] [CrossRef] [Google Scholar]
  123. White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52 [Google Scholar]
  124. White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]
  125. Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
  126. Yang, T., Cai, Y. -C., Cui, W., et al. 2022, MNRAS, 516, 4084 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Calculating baryonic component percentages

The baryonic component percentages are calculated by a fitting of a normal distribution, given by the Eq. (A.1), where the baryonic X percetage is given by %MX,bar≡(MX,bar/Mbar)×100 per cent, with MX,bar the baryonic X component mass and Mbar the total baryonic mass of the cluster. Also μ is the mean baryonic percentage value and σ os the standard deviation of it.

f ( % M X , bar ) = 1 2 π σ e − ( % M X , bar − μ ) 2 2 σ 2 $$ f({\%}M_{\mathrm {X, bar}}) = \frac {1}{\sqrt {2\pi } \sigma } e^{-\frac {\left ({\%}M_{\mathrm {X, bar}} - \mu \right )^{2}}{2 \sigma ^{2}}} $$(A.1)

thumbnail Fig. A.1.

Percentage distributions of cold-gas (T<105 K in upper row panels), warm-gas (105≤T<107 K in upper-center row panels), hot-gas (T≥107 K in bottom-center row panels) and all-gas (in bottom row panels) masses considering the baryonic masses. The histograms are divided in the haloes with log 10 ( M 500 c / M ⊙ ) < 14 $ \log _{10}{\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} < 14 $ (left colum panels), log 10 ( M 500 c / M ⊙ ) ≥ 14 $ \log _{10}{\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} \geq 14 $ (center colum panels) and all halo (right colum panels) masses. (Black histograms) TNG100 haloes. (Red histograms) TNG300 haloes. The dashed lines represent the normal distribution fitted in their respective histogram (with their respective colour).

Fig. A.1. shows the distributions of the gas component percentages (cold, warm, and hot), while Fig. A.2. illustrates the distributions of the stellar component percentages (galaxies and ICL, using rop and Subfind apertures). In this context, comparing these percentages with the total baryonic mass of the clusters is crucial, as it provides a straightforward way to analyse the scatter in TNG100 and TNG300 haloes. For instance, low-mass haloes exhibit a high scatter in the hot gas percentages, which consequently results in a high scatter across all haloes for this component. A similar trend is observed for warm and cold gas components, where the standard deviation is comparable to the mean values. This occurs due to the significant number of clusters with low percentages for these components, which affects the shape of the distribution. Nevertheless, there are numerous clusters with high warm and cold gas percentages, supporting the inclusion of these components as ”missing baryons” within the simulation environment. Additionally, the standard error, calculated as σ / N $ \sigma / \sqrt {N} $ (where N is the sample size), is on the order of tenths (see Table 3), making these percentages comparable to each other. Another baryonic component identified as missing is the ICL, which is derived from the difference between the stellar halo and the galaxies within the halo. The differences between rop and Subfind apertures are evident, with the latter capturing most of the stellar components of the haloes (see Section 4.1 and Fig. 1) for details.

thumbnail Fig. A.2.

Porcentual distributions stellar masses considering the baryonic masses. Stellar galaxy mass in rop aperture (upper row panels), ICL mass in rop aperture (upper centre row panels), stellar galaxy mass in Subfind aperture (bottom centre row panels) and ICL mass in Subfind aperture (bottom row panels). The histograms are divided in the haloes with log 10 ( M 500 c / M ⊙ ) < 14 $ \log _{10}{\left ( M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} < 14 $ (left colum panels), log 10 ( M 500 c / M ⊙ ) ≥ 14 $ \log _{10}{\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} \geq 14 $ (centre colum panels) and all halo (right colum panels) masses. (Black histograms) TNG100 haloes. (Red histograms) TNG300 haloes. The dashed lines represent the normal distribution fitted in their respective histogram (with their respective colour).

All Tables

Table 1.

Parameters of TNG100 and TNG300 simulations from IllustrisTNG (Pillepich et al. 2018c).

In the text
Table 2.

Mass notations.

In the text
Table 3.

Percentages of baryonic components.

In the text
Table 4.

Gas mass percentages.

In the text
Table 5.

Fit values of SPL, BPL and GDPL models.

In the text

All Figures

thumbnail Fig. 1.

Distribution of the percentage difference between the stellar masses from Subfind and 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2}\sim r_{\mathrm {op}} $ apertures under different constrains. The black, red, and blue lines correspond to measurements within halocentric radius R500c, R200c, and the mass considered by the FoF algorithm, respectively. The dashed lines represent the normal distribution fitted to each constraint and are shown in their corresponding colour.
In the text
thumbnail Fig. 2.

SPL fit of baryon mass fraction fbar in IllustrisTNG at R500c and redshift z. Black dots show individual haloes; the green region indicates the SPL fit. Each panel is labelled with their corresponding simulation tag and the value for the αSPL slope from the SPL fit in the upper-left corner. Top panels: The fraction of baryons includes the hot gas and the galaxy stellar mass within the 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ aperture. Bottom panels: fraction of baryons includes all the gas and all stars inside the halo. The shaded areas are represented with a width of 3σ.
In the text
thumbnail Fig. 3.

SPL fit of baryon mass fraction fbar in IllustrisTNG at R500c for all redshifts. Black dots show individual haloes; the green region indicates the SPL fit. Each panel is labelled with their corresponding simulation tag and the value for the αSPL slope from the SPL fit in the upper-left corner. Top panels: fraction of baryons includes the hot gas and the galaxy stellar mass within the 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ aperture. Bottom panels: fraction of baryons includes all the gas and all the stars inside the halo. The yellow areas represent the Planck Collaboration XIII (2016) fraction of baryons. The C18 data (blue dots) is also shown, with its SPL fit (blue areas). The shaded areas are represented with a width of 3σ.
In the text
thumbnail Fig. 4.

SPL fit (green region) of the IllustrisTNG haloes (black dots), to the scaling relation of the mass fraction of baryons fbar and halo mass M measured by the FoF algorithm, for all redshifts. Each panel is labelled with the corresponding simulation tag and the value of the αSPL slope in the upper-left corner. The baryon fraction includes all gas and stellar mass within each halo. The yellow region represents the baryon fraction from Planck Collaboration XIII (2016). All regions are in 3σ.
In the text
thumbnail Fig. 5.

SPL fit (green region), BPL fit (magenta region) and GDPL fit (red region) of the IllustrisTNG haloes (median of the data in black solid line), to the scaling relation of the mass fraction of baryons fbar and halo mass M in a radius R500c. Each panel is labelled with the corresponding simulation tag in addition to a vertical dark green dashed line that represents the Mpivot of the broken fits. The power values of αSPL SPL fits, αBPL,1 and αBPL,2 BLP fits, and αGDPL,1 and αGDPL,2 GDPL fits, with the GDPL smoothness parameter δ and Mpivot value are in the Table 5. Top panels: fraction of baryons includes the hot-gas and the 2 × R M star / 2 ∼ r op $ 2 \times R_{M_{\mathrm {star}}/2} \sim r_{\mathrm {op}} $ aperture. Bottom panels: baryon fraction includes all gas and stellar mass within each halo. The yellow region represent the Planck Collaboration XIII (2016) fraction of baryons. All the regions are in 3σ.
In the text
thumbnail Fig. A.1.

Percentage distributions of cold-gas (T<105 K in upper row panels), warm-gas (105≤T<107 K in upper-center row panels), hot-gas (T≥107 K in bottom-center row panels) and all-gas (in bottom row panels) masses considering the baryonic masses. The histograms are divided in the haloes with log 10 ( M 500 c / M ⊙ ) < 14 $ \log _{10}{\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} < 14 $ (left colum panels), log 10 ( M 500 c / M ⊙ ) ≥ 14 $ \log _{10}{\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} \geq 14 $ (center colum panels) and all halo (right colum panels) masses. (Black histograms) TNG100 haloes. (Red histograms) TNG300 haloes. The dashed lines represent the normal distribution fitted in their respective histogram (with their respective colour).
In the text
thumbnail Fig. A.2.

Porcentual distributions stellar masses considering the baryonic masses. Stellar galaxy mass in rop aperture (upper row panels), ICL mass in rop aperture (upper centre row panels), stellar galaxy mass in Subfind aperture (bottom centre row panels) and ICL mass in Subfind aperture (bottom row panels). The histograms are divided in the haloes with log 10 ( M 500 c / M ⊙ ) < 14 $ \log _{10}{\left ( M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} < 14 $ (left colum panels), log 10 ( M 500 c / M ⊙ ) ≥ 14 $ \log _{10}{\left (M_{\mathrm {500c}} / {\mathrm {M_{\odot }}} \right )} \geq 14 $ (centre colum panels) and all halo (right colum panels) masses. (Black histograms) TNG100 haloes. (Red histograms) TNG300 haloes. The dashed lines represent the normal distribution fitted in their respective histogram (with their respective colour).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.