Revisiting the Missing Baryon Problem: Baryon Redistribution Across Halo Evolution

Introduction

The cosmic baryon fraction, defined as the ratio of baryonic to total matter, is an important prediction of any cosmological model and is tightly constrained by observations of the cosmic microwave background (e.g., [74, 97]). Despite this well-established global value, its connection to the baryonic content of individual dark matter halos remains debatable. Figure 1 presents a schematic description of the problem in a rough sense. Both observations and theoretical studies have consistently shown that the baryon fraction within halos is suppressed relative to the cosmic average, particularly at galaxy and group scales, giving rise to the long-standing “missing baryon problem” (e.g., [45, 15, 8, 58, 11, 121, 79, 112, 98]).

This discrepancy arises from the complex redistribution of baryons across galaxies, the circumgalactic medium (CGM), and the intergalactic medium (IGM), which are dynamically connected. In this cycle, gas accretes from the IGM onto dark matter halos, is processed through the CGM and interstellar medium (ISM), and is subsequently redistributed by star formation and feedback-driven outflows ([132, 92, 42, 22, 114, 94, 40, 5, 134, 89]).

Supernova (SN) and active galactic nucleus (AGN) feedback processes play a central role in regulating this cycle. SN-driven winds can efficiently remove gas from low-mass halos and enrich CGM and IGM, while AGN feedback can heat and expel gas on larger scales, altering both the thermal state of the surrounding medium and the efficiency of subsequent gas accretion [30, 87, 119, 100, 60, 6, 107, 115, 81, 57] (Oppenheimer & Davé 2008).

However, observationally probing the CGM and IGM remains highly challenging due to their diffuse nature [131, 125], and direct constraints on the baryonic content of low-mass halos are still severely limited [48, 17, 18]. As a result, cosmological hydrodynamical simulations provide valuable tools for studying the baryon cycle [92].

State-of-the-art simulations such as TNG [96], EAGLE [26], SIMBA [29], and HORIZON-AGN [35] successfully reproduce many key galaxy observables, including the galaxy stellar mass function, size–mass relation, and cosmic star formation history. However, despite this success at galaxy scales, these simulations often predict significantly different baryon distributions in the CGM and IGM, reflecting the strong dependence on subgrid prescriptions for baryonic physics [114, 129, 25, 62, 104].

A widely adopted approach to characterizing these differences is to examine the baryon fraction as a function of halo mass at fixed redshift. These distributions provide a useful statistical description of baryon content and have been extensively used to compare simulations and observations [81, 5, 134, 24]. However, such population-level statistics do not unambiguously encode the evolutionary pathways of individual halos.

This limitation arises because the baryon cycle is inherently time-dependent and involves following physical processes which need to be modeled with high enough physical resolutions. Gas accretion, cooling, heating, and ejection operate on different timescales and depend sensitively on halo mass, redshift, and environment [2, 91, 81]. Since the various processes governing baryon cycling (e.g., accretion, cooling, and feedback) operate on different timescales and do not establish a quasi-steady equilibrium, interpreting the baryon fraction distribution alone is not enough to reveal the underlying physical processes.

To overcome these limitations, an approach of fine time resolution is required, in which the assembly histories of individual halos are explicitly tracked and connected to the physical processes governing baryon accretion and loss. In this work, we use the NewClusTER cosmological zoom-in simulation [55], which provides high mass, spatial, and temporal resolutions. By using merger trees and gas tracer particles, we can follow the motion of baryons across different phases and environments, enabling a detailed reconstruction of baryon fraction evolution along individual halo histories.

The primary goals of this study are twofold. First, we demonstrate that the baryon fraction distribution at fixed redshift, as usually exercised by observational studies, does not reveal the whole evolutionary pathways of halos and identify the physical origin of this discrepancy. Second, we characterize the evolution of the baryon fraction of individual halos and identify the dominant physical processes shaping its behavior across different mass regimes.

This paper is organized as follows. In § 2, we describe the simulation and methodology. In § 3, we compare the baryon fraction distribution at fixed redshift with individual evolutionary tracks and investigate their physical origin. In § 4, we explore the connection between baryon fraction scatter and large-scale environment, compare our results with observational constraints, and discuss the limitations of our analysis. Finally, in § 5, we summarize our main conclusions.

Methods

In this section, we provide a detailed description of the simulation. Since halo detection and the mechanism of tracer particles are key components of this study, we discuss them in further detail below.

Numerical simulations

We use the NewCluster simulation [55], which is a cosmological zoom-in simulation targeting a cluster-scale region. At redshift $z = 0$, the virial mass of the main galaxy cluster is $\sim10^{14.67}\,\mathrm{M}_{\odot}$, corresponding to that of a Virgo-like cluster. The simulation was carried out with RAMSES-yOMP [53], which implements an OpenMP-parallelized extension of the RAMSES code originally developed by Teyssier [123]. The simulation adopts the Chabrier (2003) initial mass function and is initialized using the WMAP7 cosmology [74]. The total comoving volume is $(142\,\mathrm{Mpc})^3$ and a zoom-in cluster region at the center with a size of $(28.5\,\mathrm{Mpc})^3$. In the zoom-in region, the best spatial resolution of gas cells is $68\,\mathrm{pc}$. The mass resolutions of stellar and dark matter particles are $2 \times10^4\,\mathrm{M}_{\odot}$ and $1.3 \times10^6\,\mathrm{M}_{\odot}$, respectively. The corresponding pseudo-lagrangian gas mass resolution is $\sim10^5\,\mathrm{M}_{\odot}$.

In NewCluster, radiative cooling and heating are implemented under the assumption of collisional ionization equilibrium. Cooling rates are evaluated using the tabulated values of Sutherland & Dopita [122], allowing primordial gas to cool down to temperatures of $10^4\,\mathrm{K}$. At lower temperatures, metal-line cooling further permits the gas to reach temperatures as low as $0.1\,\mathrm{K}$ [28]. For the reionization model, a uniform ultraviolet background is introduced at $z = 10$, following the model of Haardt & Madau [50]. Star formation is permitted in individual cells once the gas density exceeds $5\,\mathrm{H}\,\mathrm{cm}^{-3}$, following the gravoturbulent model of Kimm et al. [71]. NewCluster does not apply the turbulent Jeans length criterion for star formation, thereby preserving scale invariance as described in Dubois et al. [36]. The star formation efficiency is calculated following the equation of Federrath & Klessen [41].

Stellar feedback in the simulation is modeled through supernovae (SNe) and stellar winds. SN feedback is treated by distinguishing between type II and type Ia events (SNII and SNIa), both of which are implemented using the mechanical feedback framework developed by Kimm & Cen [70]. This scheme incorporates an empirical momentum boost that depends on the local gas density at the site of the explosion [71]. Each star particle undergoes SNII explosions after reaching an age of $5\,\mathrm{Myr}$, and the corresponding feedback energy injected per unit stellar mass is fixed to $e_{\mathrm{SNII}} = 0.045 \times10^{51}\,\mathrm{erg}\,\mathrm{M}_{\odot}^{-1}$.

Massive black holes (MBHs) are seeded in dense, evolved stellar environments based on gas and stellar properties. Seeds form in cells with gas density $n_{\mathrm{cell}} > 100\,\mathrm{H}\,\mathrm{cm}^{-3}$ and stellar mass $m_{\star,\mathrm{cell}} > 1.88 \times10^7\,\mathrm{M}_{\odot}$, with additional constraints that the stellar velocity dispersion is $\sigma_{\star} > 50\,\mathrm{km}\,\mathrm{s}^{-1}$ and stellar age $t_{\star} > 50\,\mathrm{Myr}$ to avoid short-lived clumps as the sites of black hole seeds. An exclusion radius of $4\,\mathrm{kpc}$ is adopted to prevent multiple seeds within a single galaxy, ensuring that MBHs form in stable, long-lived stellar systems.

Feedback from AGN is incorporated via a dual-mode model, distinguishing between radio/jet and quasar/heating modes depending on the accretion rate and on the spin of the MBH [34, 36]. The Eddington ratio is defined as $\dot{M}_{\mathrm{BHL}}/\dot{M}_{\mathrm{Edd}}$, where $\dot{M}_{\mathrm{BHL}}$ is the Bondi–Hoyle–Lyttleton accretion rate and $\dot{M}_{\mathrm{Edd}}$ is the Eddington accretion rate. When this ratio exceeds $0.01$, the AGN operates in the quasar mode, in which isotropic thermal energy is injected to heat the surrounding gas cells, thereby suppressing star formation in the central region. For values below $0.01$, the jet mode is activated, where momentum is imparted to the gas in the direction parallel to the spin axis of the MBH, driving outflows from the host galaxy. For theoretical consistency, the Eddington ratio is restricted to values not exceeding unity¹. Further details of the numerical implementation are provided in the overview paper by Han et al. [55].

Since NewCluster has a high stellar mass resolution of $2 \times10^4\,\mathrm{M}_{\odot}$ and a spatial resolution of $68\,\mathrm{pc}$ at $z = 0$, feedback processes can be more effectively resolved even on a cosmological scale, which is expected to provide an advantage for measuring baryon fraction. In addition, NewCluster outputs snapshots at intervals of $15\,\mathrm{Myr}$, which enables detailed tracking of halo assembly histories and baryon fraction evolution across cosmic time. This is particularly important at high redshift, where low-mass halos evolve rapidly and thus require fine temporal resolution to capture their detailed evolutionary behavior. In Figure 2, we present an outline of NewCluster at $z = 0.44$. Panel (a) represents the dark matter density map centered on the main cluster halo in NewCluster. Since NewCluster is a zoom-in simulation targeting a cluster region, its dark matter distribution exhibits filamentary structures connected to the central halo. Furthermore, since cluster regions are very dense, environmental effects are not confined to cluster subhalos but can extend to surrounding halos beyond the cluster virial radius, which may therefore exhibit properties different from those of field counterparts [44, 101, 54, 78]. Also, gas expelled by strong cluster-scale outflows can heat and redistribute the surrounding medium, potentially affecting the baryon content of nearby halos. Although NewCluster is cluster-centered, the analyzed main-halo sample extends out to $\sim10 R_{\mathrm{vir}}$ from the main cluster at $z \sim0.44$ and are therefore not restricted to the virialized cluster environment, while still sampling an overdense region.

Panel (b) shows the baryon fraction distribution in the same region as panel (a). In the main galaxy cluster region (green dashed circle), the baryon fraction is generally lower than the cosmic baryon fraction, as indicated by the blue colors, while the brightest cluster galaxy (BCG) itself exhibits a locally enhanced baryon fraction. In contrast, in the outer regions of the cluster, the baryon fraction exceeds the cosmic value (red), primarily due to gas redistribution driven by virial shocks and large-scale outflows from the main cluster [30, 9]. Panel (c) presents an example halo, where the baryon fraction exceeds the cosmic value within the galaxy, decreases in the CGM, and increases again in the IGM, reflecting the redistribution of baryons by feedback-driven ejecta from the central galaxy [15, 112, 109, 52, 108, 4, 1].

Halo (galaxy) detection and sample selection

Galaxies and halos are independently identified with the AdaptaHOP algorithm [3, 127], which locates density maxima and saddle points in the particle distribution to define preliminary halo centers and extents. The hierarchical organization of substructures is constructed using the most massive sub-node approach, in which structures connected to the most massive leaf structure are designated as the central component of a host halo, while all other embedded structures are classified as subhalos [127]. We set a minimum of 100 particles for structure identification, corresponding to mass thresholds of approximately $10^{6}\,\mathrm{M}_{\odot}$ for galaxies and $10^{8}\,\mathrm{M}_{\odot}$ for halos, respectively. With these criteria, we identify 4,452 galaxies and 64,471 halos at $z = 0.44$ (see Jeon et al. [65], for more details).

Since this study focuses on the baryon fraction of main halos to exclude orbital motion-associated environmental effects from main halos [113, 105, 23, 65], we exclude subhalos from our sample, yielding a total of 24,227 halos. The detailed analysis was performed for halos with masses larger than $2 \times10^{10}\,\mathrm{M}_{\odot}$, ensuring that each halo is resolved with at least $\sim10^{4}$ dark matter particles and thus provides numerically robust results; this selection yields a final sample of 361 halos.

Tracer particles

NewCluster adopts an adaptive mesh refinement (AMR) approach, which follows an Eulerian framework in which the gas dynamics is computed on a grid-based representation of the fluid. AMR simulations allow high spatial resolution to be dynamically concentrated in regions of interest following the gravitational collapse, and multi-scale gas dynamics across a wide dynamic range, with an efficient shock-capturing method [90, 61]. Although this approach provides computational efficiency and offers advantages in specific aspects, it has an inherent limitation: due to its grid-based, Eulerian nature, it is difficult to track the Lagrangian history of gas elements. This limitation complicates the tracking of gas flows and, in particular, makes it difficult to accurately quantify inflow and outflow masses separately as well as the gas accretion history.

To overcome these disadvantages, NewCluster adopts tracer particles based on Monte-Carlo approach which follow the gas distribution of the simulation box [46, 13]. Each tracer particle moves probabilistically in order to match the mass exchange processes in the simulation, including intercell gas fluxes, star formation, MBH accretion, and mass ejection from feedback processes. Depending on the underlying physical process, tracers can transition back and forth between gas, star, and MBH tracer states. This approach enables the reconstruction of Lagrangian histories of baryonic matter throughout the simulation (e.g., [73] (Cadiou et al. 2022)). NewCluster includes a total of 267,997,944 tracer particles, corresponding to a uniform tracer mass resolution of $5.17 \times10^{5}\ \mathrm{M}_{\odot}$ per tracer. In NewCluster, although the spatial distribution of gas tracer particles closely follows that of the gas, their finite number and relatively poor mass resolution compared to the gas (they are five times more massive than the pseudo gas mass resolution) can introduce minor discrepancies. For more detailed descriptions of the implementation and its validation, we refer the reader to Cadiou et al. [13] (see also Han et al. [55]).

Results

In this section, we compare the baryon fraction distribution with individual evolutionary tracks. Since the shape of the evolutionary tracks is a key element of this analysis, its physical origin is discussed in detail. For convenience, we hereafter denote the baryon fraction as $f_{\mathrm{b}}$.

Baryon fraction distribution vs evolutionary tracks

In this subsection, we compare the distribution of baryon fraction, $f_{\mathrm{b}}$, as a function of halo mass at a given redshift with the evolutionary tracks that trace the growth histories of massive halos selected at $z = 0.44$. To enable this comparison, it is necessary to measure $f_{\mathrm{b}}$ for each halo. In this paper, we adopt the virial radius provided by the Adaptahop output as the halo boundary for NewCluster and Horizon-AGN. We note that the results for TNG, SIMBA, and EAGLE are taken from Ayromlou et al. [5]. In our analysis, baryon fraction is measured as the total baryonic mass within the adopted halo boundary divided by the halo mass, where the halo mass is the sum of baryonic and dark matter masses.

It should be noted that NewCluster is currently at $z = 0.44$ and is still going. We use the last snapshot data at the time of writing this manuscript, while the other simulations mentioned above are analyzed at the local universe. This redshift range ($z > 0.44$) is sufficient for our purpose because it covers most of the halo mass range we are interested in, which will be clear in the following paragraphs. Moreover, as shown in Figure A1, the redshift evolution at $z < 0.44$ is small.

In Figure 3, panel (a) presents the distribution of the halo baryon fraction, for various cosmological simulations. Their key specifications are summarized in Table 1. TNG (Pillepich et al. [96]), EAGLE (Crain et al. [26]), SIMBA (Davé et al. [29]), and Horizon-AGN (Dubois et al. [35]) are full-volume simulations, providing statistically significant samples but with relatively lower resolution. In contrast, NewCluster is a zoom-in simulation, which is designed to achieve higher resolution in selected regions at the cost of a smaller sample size, which allows more detailed numerical treatments. In particular, the median effective gas cell radius in NewCluster is 2.69 kpc, compared to 15.8 kpc in TNG100, which is also a mesh-based simulation, while the mean effective cell radius of star-forming gas is 85 pc in NewCluster (as opposed to 355 pc in TNG100). This substantially-improved spatial resolution makes NewCluster particularly well suited for the present analysis.

Simulation	Box size [cMpc]	Dark matter mass resolution [$M_{\odot}$]	Stellar mass resolution [$M_{\odot}$]	Code
TNG100	110.7	$7.5 \times10^{6}$	$1.4 \times10^{6}$	AREPO
EAGLE	100	$9.7 \times10^{6}$	$1.8 \times10^{6}$	GADGET3
SIMBA	147.6	$9.6 \times10^{7}$	$1.8 \times10^{7}$	GIZMO
Horizon-AGN	142	$8 \times10^{7}$	$2 \times10^{6}$	RAMSES
NewCluster	142 (28.5 zoom-in)	$1.3 \times10^{6}$	$2 \times10^{4}$	RAMSES

In this figure, the baryon fraction of each halo is normalized by the cosmic baryon fraction of each simulation and shown as a function of halo mass. Considering the mass resolution of each simulation, halos in NewCluster are resolved with at least $10^{4}$ dark matter particles above halo masses of $10^{10}\ \mathrm{M}_{\odot}$, while the other simulations reach a similar level of numerical robustness only above $10^{11}\ \mathrm{M}_{\odot}$. Halos above these limits can therefore be regarded as numerically robust systems.

In all simulations, the baryon fraction approaches the cosmic mean toward the cluster mass scale, reflecting the increased ability of massive halos to retain hot gas within their deep gravitational potential wells. In addition, baryon fraction is also known to be sensitive to subgrid physics prescriptions, most notably feedback processes [5, 134, 24, 62]. While NewCluster and Horizon-AGN are both run with similar versions of the RAMSES code, they show notable offsets in the first peak both horizontally and vertically. These differences may reflect a combination of numerical resolution and feedback prescriptions, including the higher resolution and somewhat stronger SN feedback adopted in NewCluster.

All simulations but EAGLE show a similar overall $N$-shaped trend (rise-fall-rise). The two most recent simulations, NewCluster and TNG, exhibit similar absolute baryon fraction values and overall trends. The horizontal offset in the peak position likely reflects differences in numerical resolution and the adopted subgrid physics between the two simulations. EAGLE does not show the $N$-shaped behavior. This is likely related to the combined effect of relatively strong SN feedback and weak AGN feedback. Compared to TNG and SIMBA, the stronger SN feedback in EAGLE suppresses the baryon fraction more efficiently in low-mass halos, while the absence of a jet-mode AGN prescription [26, 110] makes AGN-driven outflows less effective in massive halos, allowing the baryon fraction to quickly recover in massive halos. The interplay between these two effects may contribute to the absence of the characteristic $N$-shaped trend [5, 134].

Overall, these comparisons suggest that variations in the baryon fraction distribution likely reflect a combination of differences in the adopted subgrid physics and numerical resolution. To further examine whether the distinct distribution in NewCluster is related to the zoom-in simulation setup, we additionally compare other RAMSES-based zoom-in simulations in Appendix A.

Given that most simulations display a characteristic $N$-shaped trend, we next investigate whether the massive halos in NewCluster sample evolve along this trend. To investigate this, we select relatively massive halos with $M_{\mathrm{halo}} > 10^{11}\,\mathrm{M}_{\odot}$ at $z = 0.44$ in NewCluster, resulting in a sample of 65 main halos. We construct merger trees using our in-house code PhantomTree (see Appendix § 7 of Jeon et al. [64] for details), identify their progenitors, and trace their evolutionary histories, from which individual evolutionary tracks are constructed.

Panel (b) of Figure 3 shows the mean evolutionary track obtained from these 65 halos (purple), overlaid with the baryon fraction distribution of the NewCluster simulation shown in panel (a). In the high-mass regime ($M_{\mathrm{halo}} > 10^{11}\,\mathrm{M}_{\odot}$), the evolutionary tracks are almost identical to the distribution. In contrast, in the low-mass regime ($M_{\mathrm{halo}} < 10^{11}\,\mathrm{M}_{\odot}$), the baryon fraction is significantly suppressed. This indicates that halos did not evolve by simply following the distribution at present-day and that their evolutionary paths differ, particularly at low masses. A key distinction between the two curves in the low-mass regime lies in their evolutionary context. Unlike the low-mass halos in the $f_{\mathrm{b}}$ distribution at $z = 0.44$, the low-mass regime of the evolutionary tracks corresponds to earlier stages through which these halos evolved, and therefore to higher redshifts.

We note that these evolutionary tracks are subject to a progenitor-selection effect, as they are constructed by selecting relatively massive halos at $z = 0.44$ and tracing their progenitors backward in time. By construction, their low-mass progenitors are therefore selected to undergo substantial subsequent mass growth and are not representative of the general low-mass halo population in terms of their mass assembly histories. However, we find no clear systematic offset in $f_{\mathrm{b}}$ between the selected progenitors and the general halo population at comparable halo mass and redshift (plot not shown). Thus, while the sample is biased in terms of mass assembly history, this selection does not appear to introduce a strong systematic bias in the $f_{\mathrm{b}}$ evolutionary tracks.

In conclusion, the difference between the evolutionary tracks and the $f_{\mathrm{b}}$ distribution at $z = 0.44$ suggests a redshift dependence in the low-mass regime, indicating that low-mass halos may experience different environments and evolutionary processes across redshift. To clearly illustrate the redshift dependence in the low-mass regime, we present Figure 4, which shows $f_{\mathrm{b}}$ as a function of halo mass at various redshifts. In the low-mass regime, baryon fraction suppression appears progressively from lower to higher halo masses as cosmic time advances. This trend indicates that baryon fraction suppression extends toward progressively more massive halos with increasing time.

Thermal state of IGM

In NewCluster, we note that the redshift dependence of baryon fraction suppression, when examined at fixed halo mass, is found to be primarily driven by a reduction in the gas mass rather than the stellar mass. This indicates a systematic decline in gas accretion efficiency at fixed halo mass toward lower redshift.

To investigate the redshift dependence of baryon fraction suppression in low-mass halos, we construct three halo samples at fixed halo mass across different redshifts. The M8, M9, and M10 samples correspond to the halos with masses of $M_{\mathrm{halo}} = (1\text{--}2) \times10^{8}\,\mathrm{M}_{\odot}$, $M_{\mathrm{halo}} = (1\text{--}2) \times10^{9}\,\mathrm{M}_{\odot}$, and $M_{\mathrm{halo}} = (1\text{--}2) \times10^{10}\,\mathrm{M}_{\odot}$, respectively, at the current $z$. The redshift dependence of the median baryon fraction is shown in panel (a) of Figure 5. The thermal state of the surrounding IGM is quantified separately in panel (b), where we present the median IGM temperature measured from gas within 1–2 $R_{\rm vir}$ of each halo. The shaded bands represent the virial temperatures of halos in each sample, estimated using the standard relation (Bryan & Norman [12]):

where $f_T\simeq0.8$ is a fudge factor for gas temperature, $\mu$ is the mean molecular weight, $\Delta_{\rm c}(z)$ is the virial overdensity with respect to the critical density of the Universe at redshift $z$, and $E^2(z)\simeq\Omega_{\rm m}(1+z)^3+\Omega_\Lambda$, where $\Omega_{\rm m}$ and $\Omega_\Lambda$ are the cosmic matter density and the cosmic dark energy density, respectively. In this work, we adopt a simplified form of this relation by assuming $\mu=0.6$, appropriate for ionized primordial gas, and since we focus on the redshift range $z\sim0.44$–10, the virial temperature reduces to

For the M8 sample, we find that baryon fraction suppression begins at $z\sim6$, whereas the M9 and M10 samples exhibit mild increases in baryon fraction persisting until $z\sim4$ and $z\sim1.6$, respectively. The baryon content is regulated by the tug of war between $T_{\rm IGM}$ and $T_{\rm vir}$ with a delayed response associated with the filtering scale. As shown in panel (b), the onset of strong suppression occurs after the IGM becomes hotter than the halo virial temperature for a sufficient period of time. It is progressively more difficult with time for the surrounding IGM gas to be accreted onto these halos beyond this point ([47, 102, 86, 85]). As lower-mass halos have shallower gravitational potential wells, their gas accretion efficiency is more sensitive to the thermal state of the surrounding medium, causing this critical transition to occur at earlier times, and therefore leading to an earlier onset of baryon fraction suppression in the M8 sample.

The close connection between baryon fraction suppression and the thermal state of the surrounding IGM suggests that IGM heating is a key driver of the reduced gas content in low-mass halos. Heating of the IGM in NewCluster arises primarily from three mechanisms: cosmic reionization, SN feedback, and AGN feedback.

(i) Reionization: since NewCluster does not explicitly include radiative transfer, reionization is implemented using a uniform ultraviolet background model. In NewCluster, reionization begins at $z=10$ and is almost completed at $z\simeq6.1$. During this process, cold gas is heated to temperatures of order $\sim10^4$ K, making it more difficult for low-mass halos to accrete it ([128, 39]).

(ii) SN feedback: in NewCluster, some halos that form early already host star formation at $z\gtrsim10$. The resulting SN explosions generate shock waves that propagate outward from the central regions, releasing energy into the surrounding medium. This process not only expels gas from halos but also heats the surrounding gas, including material located outside the halo ([31, 60]).

(iii) AGN feedback: AGN feedback is one of the most powerful feedback mechanisms operating within galaxies ([27, 38, 72]). However, in NewCluster, due to the stringent black hole seeding criteria, black holes are rarely seeded in halos with masses below $10^{11}\,\mathrm{M}_{\odot}$. As a result, AGN activity is rare at $z \gtrsim5$, and AGN feedback is unlikely to contribute significantly to the suppression of the baryon fraction at early times for our samples. Instead, AGN feedback becomes increasingly important at later epochs, particularly from $z \sim2$ onward, when many halos enter the jet-mode stage. At these redshifts, AGN-driven outflows substantially heat the surrounding IGM and are expected to suppress the baryon fraction in more massive halos.

After the onset of reionization ($z \simeq10$), the uniform ultraviolet background radiatively heats the gas to temperatures of $\sim10^{4}\,\mathrm{K}$, which can contribute to the suppression of the baryon fraction. However, as shown in Figure 5, no immediate baryon fraction suppression is observed in any of the samples.

For the M8 sample, the IGM temperature exceeds the halo virial temperature around the onset of reionization, while baryon fraction suppression emerges only near the end of reionization. This delayed suppression is consistent with the filtering scale picture and suggests that reionization-driven heating dominates the baryon fraction suppression in M8 halos. This is consistent with the radiative transfer simulations of Katz et al. [67].

Although SN feedback already exists from $z \sim10$ in early-formed galaxies, it does not increase the global IGM temperature above the $\sim10^{4}\,\mathrm{K}$ level established by reionization by the time of $z \sim6$. For the M9 and M10 samples, however, the regime where $T_{\mathrm{IGM}} > T_{\mathrm{vir}}$ is reached only after the IGM temperature has risen above the reionization heating level. Since AGN jet-mode feedback is still not dominant during the redshift interval $z \sim6$–2, the subsequent baryon fraction suppression in the M9 and M10 halos is most naturally explained by SN feedback, which also appears to be the dominant source of the gradual IGM heating during this epoch.

Toward lower redshift ($z \lesssim2$), the IGM temperature rises sharply, coinciding with the onset of widespread AGN jet-mode activity. This behavior is consistent with previous studies showing that the thermal state of the IGM at low redshift is strongly influenced by AGN feedback (e.g., Sorini et al. [115]). The associated heating of the surrounding gas is therefore expected to extend baryon fraction suppression to halos more massive than the M9 and M10 samples. Consequently, our results suggest that the dominant source of IGM heating evolves with cosmic time: reionization-driven heating primarily affects the M8 sample, SN feedback appears to dominate the gradual IGM heating during $z \sim6$–2, which is responsible for suppression in the M9 and M10 samples, and AGN feedback becomes increasingly important at lower redshift, extending suppression to higher halo masses than those affected by SN feedback alone. Our time-resolved analysis further demonstrates that this thermal evolution of the IGM is directly linked to the progressive suppression of the baryon fraction from low- to high-mass halos.

Feedback processes shaping the evolutionary tracks

In this subsection, we discuss which feedback processes primarily shape the overall trend of $f_{\mathrm{b}}$ evolutionary tracks. Figure 6 presents the mean $f_{\mathrm{b}}$ evolutionary track of the halo sample (with $M_{\mathrm{halo}} > 10^{11}\,\mathrm{M}_{\odot}$ at $z = 0.44$; purple line), together with the evolutionary track of a selected target halo (black line). The target halo broadly follows the mean trend, indicating that it is representative of the typical evolutionary behavior in the sample.

The gray dashed lines mark the two characteristic local extrema along the evolutionary tracks. After the peak, the baryon fraction decreases, whereas beyond the trough it increases again. Based on these features, the evolutionary history of halos can be divided into three distinct phases.

First phase: cold gas accretion

During the first phase, corresponding to the increasing trend of the baryon fraction, dark matter collapses first, establishing the gravitational potential well, followed by the accretion of baryons into this pre-existing potential ([133]). The nearly monotonic increase in the baryon fraction reflects the fact that the halos in our sample form mainly at high redshift, when feedback processes are still weak or even absent. At these epochs, the gas is predominantly cold, allowing efficient baryon accretion with little thermal resistance.

Second phase: SN-domination

Once the halos grow beyond a mass of $M_{\rm halo} \sim10^{10}\,M_\odot$, corresponding to the peak in the evolutionary tracks, they enter a second phase in which the baryon fraction begins to decline. This decrease can be understood as the combined outcome of two processes: (i) depletion of baryons from halos driven by strong outflows, and (ii) inefficient baryon accretion relative to dark matter growth, leading to a reduced baryon fraction.

Several physical mechanisms can, in principle, inhibit gas accretion, including mergers, virial shocks, SN feedback, and AGN feedback. However, the halos in this second phase are not massive enough to develop stable virial shocks ([30]). Moreover, in NewCluster, due to the adopted black hole seeding criteria, black holes are not typically seeded in halos with $M_{\rm halo} \lesssim10^{11}\,M_\odot$. As a result, AGN feedback is largely absent in halos below this mass threshold. Even after a black hole is formed, AGN activity remains dominated by the quasar mode in halos with $M_{\rm halo} \lesssim10^{12}\,M_\odot$ in NewCluster where feedback primarily acts as isothermal heating that suppresses central star formation ([124]) rather than driving large-scale gas outflows (jet mode). Consequently, AGN feedback is unlikely to play a major role in shaping the baryon fraction evolution during this phase.

This interpretation is further supported by the NewHorizon2 simulation (Yi et al. 2024), which adopts almost identical subgrid feedback physics to NewCluster but suffers from unstable black hole settling problem that prevents many galaxies from retaining central black holes (see Figure A2). Despite this, the $f_{\rm b}$ distributions for halos with $M_{\rm halo} \sim10^{10}$–$10^{12}\,M_\odot$ are remarkably similar between NewHorizon2 and NewCluster (see Figure A1 in Appendix § 7). This close agreement supports the interpretation that the characteristic shape of the evolutionary tracks during the second phase is not primarily driven by AGN feedback, leaving SN feedback as the most plausible mechanism within the adopted NewCluster model.

We utilize the tracer particles in order to investigate the role of outflows in $f_{\rm b}$ depletion, overcoming the limitation of the Eulerian approach. We select tracer particles associated with the target halo which is the same as the one in Figure 6, based on their minimum radial distances. In Figure 7, panel (a) shows the spatial distribution of tracer particles at the end of the second phase, selected as those that resided within the halo virial radius at least once while the halo mass grew from $10^{10}$ to $10^{12}\,\mathrm{M}_{\odot}$. A significant fraction of these tracers lie outside the virial radius (red circle), indicating significant influence of outflows. Quantitatively, approximately 18% of the tracer particles escape beyond $R_{\mathrm{vir}}$ after reaching $r < R_{\mathrm{vir}}$.

However, as shown by the evolutionary track of the target halo in Figure 6, the baryon fraction decreases from a value close to the cosmic $f_{\mathrm{b}}$ at the beginning of the second phase to roughly half of that value by its end. The removal of only 18% of baryons by outflows is therefore insufficient to fully account for the observed decline in the baryon fraction, indicating that outflows alone cannot fully explain the suppression.

To capture additional processes that contribute to the reduction of baryon fraction, we present panel (b). This panel shows tracer particles that satisfy two conditions: they have entered within $2R_{\mathrm{vir}}$ at least once during the second phase, but have never crossed within $R_{\mathrm{vir}}$. These tracers therefore represent gas residing exclusively in the halo’s surrounding environment. Their spatial distribution at the end of the second phase reveals that a non-negligible fraction of such surrounding gas is also displaced outward, with approximately 13% of these tracers found beyond $2R_{\mathrm{vir}}$.

To investigate the connection between these baryon-flow changes and star formation activity, as well as the role of preventative feedback ([80, 81, 134]), we analyze the coupled evolution of baryon fraction, star formation rate, and baryon inflow/outflow in the target halo (Figure 8). Although the average evolutionary tracks in the second phase appear to decline monotonically, the individual target halo in fact experiences a more episodic evolution. While the overall trend is decreasing, we identify three distinct intervals during which the baryon fraction drops rapidly; these intervals are highlighted as gray shaded regions.

In panel (a), the black line shows the evolutionary track of the selected target halo, previously introduced in Figure 6, while the yellow line represents the star formation rate of the central galaxy averaged over the past 30 Myr. During the gray shaded intervals, the baryon fraction decreases while the star formation rate simultaneously increases, suggesting a close connection between baryon fraction suppression and star formation-driven feedback. In the first shaded region, a major merger event occurs. This merger triggers star formation, leading to baryon loss through feedback ([118, 59]). In addition, while

dark matter is relatively unaffected by stellar feedback and continues to accrete efficiently during the merger, baryons experience significant resistance ([92]). As a result, the baryon fraction drops sharply, followed by a partial recovery as baryons accrete at later times. This behavior suggests that mergers can induce rapid, transient variations in the baryon fraction, but do not significantly alter the overall declining trend in the second phase.

Panel (b) presents the inflow (blue) and outflow (red) masses, inferred from tracer particles, measured at $R_{\mathrm{vir}}$ in the corresponding mass bins along the evolutionary track. During the gray shaded intervals, strong outflows are clearly observed, coinciding with the enhanced star formation activity seen in panel (a). At the same time, inflow is suppressed relative to the baryonic mass expected from halo growth under the assumption of the cosmic baryon fraction, $\Delta M_{\mathrm{halo}} \times f_{\mathrm{b,c cosmic}}$, during these periods. The combination of enhanced outflows and reduced inflows is consistent with a preventative-feedback scenario, in which feedback not only removes gas from within the halo but also inhibits the accretion of surrounding gas. As also shown in panel (b) of Figure 7, this process can even expel gas residing in the surrounding environment.

Taken together, the tracer results directly demonstrate enhanced outflows and reduced inflows during the declining phase. The temporal association of these changes with enhanced star formation, together with the limited contribution of AGN feedback in this mass regime under the adopted black hole model, suggests that SN feedback is likely the dominant mechanism within NewCluster for expelling gas from halos and suppressing gas accretion. Our interpretation that SN feedback dominates in this mass range is consistent with previous studies based on statistical distributions at fixed redshift [84, 56, 5, 24], while our analysis provides a more direct evolutionary perspective by tracing halo assembly histories. We note, however, that the quantitative strength of the resulting baryon-fraction decrease likely depends on the adopted subgrid physics.

Third phase: hot gas accretion

When the halo mass exceeds $\sim10^{12}\,\mathrm{M}_{\odot}$, a third phase begins, in which $f_{\mathrm{b}}$ gradually recovers. This recovery implies an increase in the gas content relative to the dark matter mass, indicating a substantially improved gas retention capability in massive halos [69, 30].

During the preceding second phase, SN-driven feedback expels gas from halos and suppresses subsequent inflows. The affected gas is predominantly heated to high temperatures, making it difficult to remain gravitationally bound and thus causing it to reside beyond the virial radius. This behavior is evident in Figure 2, panels (b) and (c), where the baryon fraction within the CGM ($R_{\mathrm{gal}} \lesssim r \lesssim R_{\mathrm{vir}}$) is lower than the cosmic value, while the surrounding IGM ($> R_{\mathrm{vir}}$) exhibits an excess of baryons. Such redistribution of baryons provides a natural explanation for the missing baryon problem observed in halos [11].

Once halos grow to masses of $\gtrsim10^{12}\,\mathrm{M}_{\odot}$, however, they become capable of sustaining hot-mode accretion. The deepened gravitational potential enables the accretion and retention of hot gas that was previously unable to cool or bind to the halo and thus remained in the IGM [77, 93, 37]. As a result, baryons that were heated and displaced by earlier feedback processes can be reincorporated into the halo, leading to the observed recovery of the baryon fraction during the third phase. This trend, in which increasingly massive halos approach the cosmic baryon fraction as their potential wells deepen toward group and cluster scales, is consistent with both previous numerical studies and observational constraints [48, 17, 18, 5, 24].

To investigate the evolution underlying this picture, we compare the radial density profiles of dark matter, gas, and hot gas ($T \ge10^{6}\,\mathrm{K}$) for the target halo shown in Figure 6 at two representative epochs, as summarized in Figure 9. Panel (a) compares the radial density profiles of dark matter, gas, and hot gas for the target halo at two epochs: at the end of the second phase, when the halo mass reaches $M_{\mathrm{halo}} \simeq10^{12}\,\mathrm{M}_{\odot}$ ($z = 3.6$, dashed lines), and at a later time when the baryon fraction has been substantially recovered, $M_{\mathrm{halo}} \simeq5 \times10^{12}\,\mathrm{M}_{\odot}$ ($z = 1.6$, solid lines).

The dark matter density profiles show little difference between the two stages when expressed in units of $R_{\mathrm{vir}}$, indicating that the underlying gravitational potential structure remains largely unchanged. In contrast, the gas density profile exhibits a noticeable enhancement in the circumgalactic region, approximately $0.25 \lesssim r/R_{\mathrm{vir}} \lesssim1$, implying that the increase in the baryon fraction is primarily driven by the accumulation of gas in the CGM.

Moreover, the hot gas component shows a significant increase across nearly the entire radial range, suggesting that the halo becomes increasingly efficient at retaining and accreting hot gas. This supports the interpretation that the recovery of the baryon fraction in the third phase is closely linked to the onset of efficient hot-mode accretion in massive halos.

However, the increase in hot gas within the CGM does not by itself prove that hot-mode accretion is the only source. It could also arise from feedback processes originating in the central galaxy. Disentangling these two effects is therefore essential. We decompose the hot gas profile at $z = 1.6$ according to its origin and address in panel (b). The hot gas is classified into three categories: (i) gas originating from the central region (“from inside”), (ii) accreted gas that was already hot at the time of infall ($T \ge10^{6}\,\mathrm{K}$), and (iii) accreted gas that was cold at infall ($T < 10^{6}\,\mathrm{K}$). Here, inflowing gas is defined using tracer particles that entered the halo after $z = 3.6$. Tracers that were already present at $z = 3.6$ and do not belong to category (i) are not included in this classification; however, they constitute only about 2% of the total hot gas and are therefore negligible.

To identify hot gas expelled by internal feedback, the galactic region is defined as the inner $0.1\,R_{\mathrm{vir}}$ (e.g., Kocjan et al. 2024). Tracer particles that passed through this region at least once are classified as “from inside.” For the remaining tracers, the gas temperature at the moment they crossed $R_{\mathrm{vir}}$ is used to distinguish between hot and cold infall. Hot gas that was already hot at infall is interpreted as evidence for the halo’s ability to accrete hot gas owing to its deepened potential well, whereas hot gas that entered cold and subsequently heated, as well as gas originating from the inner region, is attributed to feedback-related processes and virial heating in sufficiently massive halos.

As shown in panel (b), the increase in hot gas density in the radial range $0.25 \lesssim r/R_{\mathrm{vir}} \lesssim1$ is dominated by gas that was already hot at infall. This provides direct evidence that the enhanced hot gas content in the CGM originates predominantly from hot-mode accretion, rather than from in-situ heating by feedback processes.

In this mass regime ($M_{\mathrm{halo}} > 10^{12}\,M_{\odot}$), most galaxies host central black holes with properties broadly consistent with observational constraints ([55]). Also, AGN jet mode feedback is activated in most halos of this mass range and redshift ([33, 130, 51]). Since the onset of the jet mode depends primarily on the Eddington ratio rather than halo mass itself, the characteristic halo mass at which the jet mode is triggered varies from halo to halo. Once activated, the jet mode can drive strong

outflows, leading to a temporary reduction in the halo baryon fraction. In fact, such a decrease was visible in Figure 6 at $M_{\mathrm{halo}} \gtrsim5 \times10^{12}\,M_{\odot}$. However, because the activation mass and timing differ among individual halos, these episodic decreases are largely averaged out in the evolutionary tracks and do not appear as a coherent trend. Nevertheless, without sufficiently effective AGN feedback, the increase of the baryon fraction during the third phase becomes steeper, as seen for NewHorizon2 in Figure A1.

Therefore, in the third phase, the recovery of the baryon fraction is primarily driven by the halo’s increased mass and deepened potential well, which enable the accretion of gas that was already hot at infall. Our results show that this externally supplied hot gas dominates the CGM, establishing hot-mode accretion as the primary origin of the recovered baryon content. While AGN jet-mode feedback can induce decreases in the baryon fraction, these effects are transient and sub-dominant.

Discussion

Comparison with observations on galaxy scale

As direct observational constraints on halo-scale baryon fractions remain limited, especially in the low-mass regime, we instead compare galaxy-scale properties of our simulation with observations, where measurements are more robust. While agreement at the galaxy scale does not guarantee consistency at the halo scale, it provides a meaningful observational baseline for assessing the credibility of the baryonic physics realized in the simulation in the absence of direct halo-scale constraints.

To this end, we compare our NewCluster sample with observational measurements from the MAGPI (Middle Ages Galaxy Properties with Integral Field Spectroscopy) survey (Foster et al. [43]). MAGPI is an integral-field spectroscopic survey conducted with MUSE on the VLT, providing spatially resolved kinematics of stars and ionized gas for galaxies at intermediate redshift, with a median redshift of $\langle z\rangle\sim0.35$, corresponding to look-back times of $\sim3$–$4$ Gyr.

Using the MAGPI data, Sharma et al. [111] constructed a sample of star-forming galaxies with stellar masses of $10^{8}$–$10^{11.5}\,M_{\odot}$, comprising 397 galaxies over the redshift range $z = 0.1$–$0.85$, and measured the dark matter fraction on galaxy scales.

Dark matter fraction $f_{\mathrm{DM}}$ within the effective radius was estimated by decomposing the dynamical mass as

where $M_{\mathrm{disk}} + M_{\mathrm{bulge}}$ were derived from spectral energy distribution fitting, $M_{\mathrm{HI}}$ and $M_{\mathrm{H}_{2}}$ gas masses were inferred using scaling relations, and $M_{\mathrm{dyn}}$ was obtained from rotation curve fitting. The factor of 1.33 in HI gas accounts for the contribution of helium to the atomic gas mass. For further methodological details, we refer the reader to Sharma et al. [111].

For a fair comparison with MAGPI, we construct a galaxy sample from the last snapshot of NewCluster by applying selection criteria as closely matched as possible to the observations. The redshift of the last snapshot of NewCluster (0.44) is reasonably close to the median redshift of the MAGPI data. Since satellite galaxies are strongly affected by environmental processes, we restrict our analysis to central galaxies only. We further apply the same stellar mass cut as in the MAGPI sample, selecting galaxies with $10^{8} \le M_{\star}/M_{\odot} \le10^{11.5}$. In addition, the MAGPI analysis accounts only for cold gas when estimating the baryonic mass. To ensure consistency, we therefore consider only gas with temperatures below $10^{4}\,\mathrm{K}$, corresponding to the HI temperature regime, and densities above $1\,\mathrm{H\,cm^{-3}}$ in the NewCluster sample. Based on these criteria, we measure the stellar mass and cold gas mass within the half-mass radius and define the baryonic mass as

In the approach adopted by Sharma et al. [111], the dark matter mass is inferred by subtracting the stellar and cold gas components from the dynamical mass, such

that warm or hot gas components are not explicitly removed, even though they are implicitly included in the dynamical mass estimate. To apply the same limitation consistently to NewCluster, we compute this pseudo dark matter mass by subtracting $M_{\mathrm{baryon}}$ from the total matter mass within the half-mass radius.

Figure 10 presents a comparison between MAGPI and NewCluster. The $x$-axis shows the baryonic surface density within the effective radius (or the half-mass radius for NewCluster), while the $y$-axis indicates $f_{\mathrm{DM}}$ within the corresponding region. We find that the NewCluster galaxies broadly match the observed galaxy-scale $f_{\mathrm{DM}}$ measured in the MAGPI sample. Although $f_{\mathrm{b}}$ in the CGM is known to be highly sensitive to subgrid feedback prescriptions, galaxy-scale baryon and dark matter fractions are also directly affected by these processes. The agreement with MAGPI therefore provides an important observational consistency check on the galaxy-scale baryonic physics realized in NewCluster and a meaningful basis for investigating the baryon cycle on larger scales. We emphasize, however, that agreement on galaxy scales does not by itself guarantee an accurate baryon distribution at the halo scale.

Environmental effects

While the overall shape of the baryon fraction evolutionary tracks appears to be primarily governed by internal processes and surrounding gas conditions, a significant halo-to-halo scatter is present at fixed halo mass. In particular, halos with similar masses can exhibit substantially different baryon fractions during their evolution. Musso et al. [83] used excursion theory to explain how the geometry and connectivity of the cosmic web influence halo assembly and the statistics of matter accretion. Their study may provide the origin for the large scatter in $f_{\mathrm{b}}$ we found.

Here, we focus on $z = 0.8$, corresponding to the epoch just before cluster-scale mergers become dynamically important in NewCluster. This choice allows us to minimize the impact of merger-driven interactions and better isolate the role of the filamentary cosmic web environment in shaping the scatter of baryon fraction.

To quantify this diversity, we classify halos according to their deviation from the mean $f_{\mathrm{b}}$ at a fixed mass. Figure 11 shows the distribution of $f_{\mathrm{b}}$, where halos with baryon fractions greater than $1\sigma$ above the mean are defined as “baryon-rich” (hereafter BR; 62 halos), those within $1\sigma$ as inside $1\sigma$, and those more than $1\sigma$ below the mean as “baryon-poor” (hereafter BP; 66 halos). To ensure sufficient numerical resolution, we restrict our analysis to halos containing more than $\sim10^{4}$ dark matter particles, corresponding to a mass cut of $M_{\mathrm{halo}} \gtrsim2 \times10^{10}\,\mathrm{M}_{\odot}$. Halos with $M_{\mathrm{halo}} \gtrsim2 \times10^{12}\,\mathrm{M}_{\odot}$ are excluded due to their limited sample size, which prevents a robust separation into BR and BP.

Even within a narrow halo mass range, the BR and BP samples can differ in baryonic mass by more than a factor of two. Such a large dispersion suggests that in addition to internal processes and surrounding gas, more external factors may also contribute to shaping the baryon content of halos.

Since our sample consists exclusively of main halos rather than satellites embedded in groups or clusters, strong local environmental processes such as ram-pressure stripping (e.g., [68, 20, 66, 63, 106]) are unlikely to be the dominant driver of the observed scatter. It is therefore more reasonable to investigate the influence of the large-scale cosmic environment, particularly the connection to filamentary structures of the cosmic web (Tuominen et al. [126]). Halos located closer to filaments may experience different gas supply, thermal conditions, and dynamical interactions (e.g., [10, 69, 95, 21, 75]).

To characterize the cosmic web structure, we employ the publicly available Discrete Persistent Structure Extractor (DisPerSE; [116, 117]). DisPerSE is a topological structure-finding algorithm designed to identify multiscale features in a density field. Indeed, Musso et al. [83] showed that topology is the correct metric for understanding halo growth. They showed how the connectivity and local curvature of the cosmic web (as treated in an extended excursion theory framework) dictate the ‘Lagrangian history’ of galactic assembly. This approach enables a robust identification of filamentary structures across a wide range of scales.

To extract the filamentary structure, we first construct a density field within the central region of NewCluster, corresponding to a comoving volume of $(42.6\,\mathrm{cMpc})^{3}$. The density field is generated on a uniform grid with a resolution of approximately $166$ kpc using the cloud-in-cell (CIC) interpolation scheme. To reduce small-scale noise and obtain a smoother representation of the large-scale structure, we further apply an additional smoothing procedure to the gridded density field.

DisPerSE is then applied to the resulting smoothed density field. We adopt a cut value of $0.7$, which corresponds approximately to a persistence ratio equivalent to a $4\sigma$ significance level. This threshold effectively suppresses spurious structures arising from sampling noise while retaining the dominant filamentary network, and corresponds to a regime where the large-scale structural distribution becomes well separated from noise-induced features.

In Figure 12, panel (a) shows the filamentary structure identified by DisPerSE. Gray dots represent the spatial distribution of all halos, while red and blue symbols indicate the BR and BP halos, respectively. For each halo, we measure the shortest perpendicular distance to the nearest filament, and the resulting distance distributions are presented in panel (b). In this analysis, we exclude halos located within $3 R_{\mathrm{vir}}$ (4.26 cMpc) of the main cluster, since the environmental effects of the massive cluster cannot be neglected in such regions. After this selection, 52 BR and 57 BP halos remain.

In panel (b), it is clear that BR halos are found closer to filaments, compared to BP halos. These results demonstrate a correlation between halo baryon fraction and proximity to filamentary structures. A two-sample Kolmogorov–Smirnov (KS) test performed on the $D_{\mathrm{fil}}$ distributions in panel (b) confirms that the BR and BP halo populations are statistically distinct at high significance, yielding $D_{\mathrm{KS} } = 0.525$ and $p = 2.1 \times10^{-7}$.

However, one important limitation of the DisPerSE-based filament identification is that the extracted filament structure does not provide information on filament width. As a result, halos located at the same perpendicular distance from a filament may reside in physically different environments depending on the filament density structure. To account for this effect, we additionally measured cylindrical radial dark matter density profiles around individual filament segments and estimated the overdensity associated with each halo based on its radial position relative to nearby filaments. We note that the resulting overdensity trends are broadly consistent with those inferred from filament distance alone.

These results indicate a correlation between halo baryon fraction and the global filamentary environment. We note, however, that proximity to filaments does not necessarily imply more efficient gas accretion, since both gas density and temperature increase toward filament spines. In the mass range considered here ($M_{\mathrm{halo}} \gtrsim2 \times10^{10} M_{\odot}$), where redshift-dependent baryon-fraction suppression is very weak, BR halos reside closer to filaments despite their relatively hot surroundings. In an additional test (plot not shown) extending to lower halo masses, where suppression becomes stronger, the BR–BP difference in filament proximity weakens, while BP halos tend to reside in hotter IGM environments, suggesting an increasing importance of the thermal state of the surrounding gas.

We therefore interpret the present result as an environmental correlation rather than direct evidence that filament proximity enhances gas accretion. Establishing the physical connection between filamentary environment and baryon acquisition requires a detailed analysis of individual gas accretion histories, which is beyond the scope of the present study and will be addressed in future work.

Caveats

One of the key uncertainties in this work arises from the adopted subgrid physics, which plays a central role in regulating the baryon cycle. As discussed in the introduction paper of NewCluster [55], the simulation exhibits a tendency to overestimate stellar masses in low-mass halos, particularly around $M_{\mathrm{halo}} \sim10^{11}\ M_{\odot}$, as indicated by the stellar-to-halo mass relation [82, 76, 7]. Elevated stellar masses can lead to stronger SN feedback, which may in turn affect the predicted halo baryon fraction in this mass regime. This suggests that the quantitative strength of the baryon-fraction decrease during the second phase likely depends on the adopted feedback model.

In addition, the black hole seeding model in NewCluster adopts stringent criteria, resulting in black holes being rarely seeded in halos with $M_{\mathrm{halo}} \lesssim10^{11}\ M_{\odot}$. Consequently, the impact of black hole feedback cannot be assessed in this lower-mass regime. However, some of these limitations are not expected to qualitatively alter our main conclusions. Although black holes are not efficiently seeded below $10^{11}\ M_{\odot}$, feedback from AGN is unlikely to play an important role in regulating the baryon fraction in such low-mass halos.

In addition, since NewCluster is a zoom-in simulation, the statistical sample size is inherently smaller than those of typical full-volume cosmological simulations, which may lead to reduced statistical robustness and render cosmic variance non-negligible. Moreover, because NewCluster follows the formation of a massive cluster environment, it may not fully represent the average environment of the Universe. Despite that, as presented in Appendix § 7, the difference between NewCluster and NewHorizon2 in terms of the baryon fraction trend is not large. Our additional comparison with the full-volume Horizon-AGN simulation, also presented in Appendix § 7, shows that low-mass halos in cluster environments have somewhat lower baryon fractions than those in non-cluster environments, while the overall trend remains similar.

Conclusions

In this work, we investigate the evolution of the baryon fraction within dark matter halos using the NewCluster cosmological zoom-in simulation, with the aim of linking baryon fraction evolution to the underlying physical processes governing halo growth. Rather than simply connecting evolutionary pathways solely from comparisons of $f_{\mathrm{b}}$ distributions at different redshifts, we directly trace individual halo histories using merger trees and tracer particles, thereby connecting baryon fraction evolution to the baryon cycle.

The high mass and spatial resolution of NewCluster, combined with its short snapshot intervals, enables a detailed investigation across cosmic time, while the use of tracer particles allows us to track gas motions and origins, making NewCluster particularly well suited for both statistical characterization of the $f_{\mathrm{b}}$ distribution and reconstruction of individual evolutionary tracks.

We trace the evolutionary histories of 65 main halos with $M_{\mathrm{halo}} > 10^{11}\ M_{\odot}$ at $z = 0.44$. Our analysis reveals an $N$-shaped trend in which the baryon fraction evolution proceeds through three distinct phases.

1. First phase: it corresponds to the early cosmic times, during which efficient cold accretion allows relatively cold gas to be readily incorporated into assembling dark matter halos.

2. Second phase: the baryon fraction appears to decline primarily due to SN feedback. Tracer analysis directly reveals enhanced outflows and reduced subsequent inflow, indicating a significant role of preventative feedback, while their association with enhanced star formation supports SN feedback as the dominant mechanism in this regime.

3. Third phase: in more massive halos, we identify a distinct recovery phase, in which deeper gravitational potentials enable efficient gas retention and reaccretion, driving the baryon fraction back toward the cosmic value. Tracer analysis reveals that this recovery is dominated by gas that was already hot at infall, demonstrating that hot-mode accretion plays a key role in restoring the baryon content. This highlights a clear transition in the dominant baryon cycle and represents a key result of this work.

We find that the $f_{\mathrm{b}}$ distribution at fixed redshift does not unambiguously encode the evolutionary pathways of halos. Halos with similar masses and baryon fractions at a given epoch can follow markedly different assembly histories, particularly in the low-mass regime where strong redshift dependence is observed. This suggests that population-level statistics alone are insufficient to capture the time-dependent nature of the baryon cycle. In particular, a direct comparison between the $f_{\mathrm{b}}$ distribution and individual evolutionary tracks reveals a systematic shift in the location of the first peak, implying that the distribution alone may misrepresent the underlying physical trends.

We further find that baryon fraction suppression in low-mass halos appears to be linked to the thermal state of the surrounding IGM. In particular, strong suppression emerges after the IGM temperature becomes higher than the halo virial temperature for an extended period of time, highlighting the critical role of environmental heating in regulating baryon accretion. In addition, the scatter in baryon fraction at fixed halo mass is correlated with the large-scale environment, with baryon-rich halos tending to reside closer to filamentary structures in mass regimes where thermal suppression remains weak.

Overall, our results demonstrate that baryon fraction evolution is governed by a complex interplay between internal feedback processes and environmental effects, and cannot be fully understood without sufficient spatial resolution and a time-resolved perspective on halo growth. These findings have important implications for the interpretation of baryon fraction measurements and for comparisons between simulations. Given that the halo baryon fraction is very sensitive to the adopted subgrid models, continued improvements and systematic comparisons across different simulations are essential to robustly uncover the physical origin of baryon distributions in the Universe.

Acknowledgements

We thank the referee for providing comments that helped clarify several points in the early manuscript. Special thanks go to Gauri Sharma for providing the MAGPI observational data in Figure 10. This work was granted access to the HPC resources of KISTI under the allocations KSC-2021-CRE-0486, KSC-2022-CRE-0088, KSC-2022-CRE-0344, KSC-2022-CRE-0409, KSC-2023-CRE-0343, KSC-2024-CHA-0009, and KSC-2025-CRE-0031 and of GENCI under the allocations A0150414625 and A0180416216. The large data transfer was supported by KREONET, which is managed and operated by KISTI. S.K.Y. acknowledges support from the Korean National Research Foundation (RS-2025-00514475 and RS-2022-NR070872). J.R. acknowledges support from the Institut de Physique des deux infinis of Sorbonne Université and by the ANR grant ANR-19-CE31-0017 of the French Agence Nationale de la Recherche. E.C. and C.S. acknowledge support from the Korean National Research Foundation (RS-2023-00241934). S.J. acknowledges support from the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2026-1-831-03) supervised by the Ministry of Science and ICT.

Baryon fraction distributions in RAMSES simulations

In this appendix, we compare the $f_b$ distributions of several RAMSES-based simulations. Figure A1 shows the $f_b$ distributions of NewCluster (blue), NewHorizon2 (orange), and YZiCS (black) at $z = 0.44$ (solid) and $z \sim0$ (dashed).

For YZiCS [19], the distribution becomes similar to that of SIMBA, suggesting that the distinct distribution in NewCluster is not simply a consequence of the zoom-in setup. Both YZiCS and NewHorizon2 show little evolution in $f_b$ between $z = 0.44$ and $z = 0$, indicating that the NewCluster distribution at $z = 0.44$ is broadly representative of the low-redshift regime.

To further assess the possible impact of the overdense environment sampled by NewCluster, we additionally use the full-volume Horizon-AGN simulation to compare halos in cluster and non-cluster environments at $z \sim0$, defining the former as regions within 14 Mpc of 12 halos with $M_{\rm vir} > 10^{14}\,M_\odot$, approximately matching the characteristic radial extent of the NewCluster zoom-in region. As shown in Figure A1, low-mass halos in cluster environments tend to have somewhat lower baryon fractions than their counterparts in non-cluster environments up to $M_{\rm halo} \sim10^{11}M_\odot$, while the overall trend remains similar.

This environmental suppression is qualitatively consistent with the slightly lower baryon fractions of NewCluster relative to NewHorizon2 at low halo masses. However, the magnitude of the environmental suppression is very weak, indicating that the large-scale environment is unlikely to be the primary origin of the overall trend.

The offset in the peak between YZiCS and NewCluster at $z = 0.44$ is likely driven not only by differences in subgrid physics, but also by differences in the dark matter mass resolution. In the case of NewHorizon2, the distribution at $z = 0.44$ closely follows that of NewCluster. In particular, the similarity in the mass range $10^{10}$–$10^{12}\,M_\odot$, despite the absence of well-centered black holes and their systematically undermassive nature (see Figure A2), indicates that AGN feedback is not the dominant process in this regime. Also, since black holes in NewHorizon2 are not well centered within their host halos, AGN feedback cannot operate efficiently. This likely causes the third evolutionary phase to begin at slightly lower halo masses and leads to weaker suppression of gas inflow, resulting in a steeper increase in baryon fraction with halo mass.