Modeling Galaxy Formation in Cosmological Simulations with CRK-HACC

TL;DR

This work extends the CRK-HACC framework with a calibrated suite of subgrid models for radiative cooling, star formation, chemical enrichment, winds, and AGN feedback to enable survey-scale, self-consistent hydrodynamic cosmological simulations. By leveraging GPU acceleration and exact integration schemes, the authors achieve exascale-ready performance while maintaining numerical stability and reproducibility across parallel domains. Calibrated against the galaxy stellar mass function (0 < z < 2) and massive-cluster gas-density profiles, the fiducial model produces realistic global star formation histories, galaxy sizes, metallicities, black hole growth, and halo gas fractions, with quantified resolution limitations in low-mass, high-z regimes. The framework lays the groundwork for large-volume mock surveys and cross-probe analyses, while highlighting the ongoing need to refine subgrid physics and resolution-variance effects as computational capabilities grow.”

Abstract

Self-consistently modeling baryonic effects in survey-scale cosmological simulations has become increasingly important as the diversity, precision, and statistical reach of modern observations continue to improve. The advent of exascale computing now enables a new generation of simulations that couple these physical processes across full-sky volumes with excellent statistical sampling of large-scale structure tracers such as galaxies, groups, and clusters. To support these efforts, we extend the CRK-HACC framework, a GPU-accelerated cosmological hydrodynamics code, with a suite of astrophysical subgrid models that simulate radiative cooling, star formation, stellar evolution, and AGN feedback within a numerically robust formulation optimized for scalability on modern exascale architectures. The models were selected and calibrated to reproduce observed galaxy stellar mass functions over the redshift range $0 < z < 2$ and cluster populations probed by cosmological surveys, capturing the large-scale baryonic evolution relevant for multi-wavelength, cross-correlated analyses. We describe the implementation and calibration of these models and demonstrate their consistency with observed galaxy population statistics and modern hydrodynamic simulations, establishing the baseline for exascale efforts that extend this framework to survey-scale volumes.

Figures (16)

• Figure 1: Cumulative distribution function of the relative difference in cooling rates computed with CLOUDY when using individual versus solar-scaled metal abundances. The individual-abundance cooling rates are calculated for one million randomly selected gas particles from the m12a FIRE-3 simulation hopkins2023firesultan2025 that tracks nine metal elements. These are compared to rates derived assuming solar-scaled abundances. Square brackets in the legend indicate the metallicity bin of each curve with the round brackets denoting the fraction of CRK-HACC gas particles in that bin. Vertical dashed lines show the median error of each bin. Among metal-enriched gas particles with $Z/Z_\odot \geq 10^{-4}$, we find that $93\%$ ($58\%$) have relative errors below $10\%$ ($1\%$). This translates into $98\%$ ($90\%$) of all simulation particles, confirming that nearly all gas in CRK-HACC cools at rates consistent with element-by-element tracking.
• Figure 2: Cumulative mass lost for the total (top), helium (middle), and metal (bottom) mass components of an SSP with solar metallicity, a primordial helium fraction, and initial mass $M_*$. Blue and red curves separate the mass-loss contributions from supernovae and stellar winds, respectively, and the black curve shows their sum. Solid lines show results from the integrated method using a timestep of $\unit{10}{\Myr}$, representative of the largest hydrodynamic steps in CRK-HACC. Shaded lines indicate the baseline solution obtained by summing instantaneous FIRE-3 rates over very fine timesteps, while dashed lines show the discretized case, where instantaneous rates are sampled only at $\unit{10}{\Myr}$ intervals. Vertical dotted blue lines mark 3.7 and $\unit{44}{\Myr}$, the characteristic timescales for core-collapse and Type Ia supernovae, and vertical dotted red lines mark 20 and $\unit{800}{\Myr}$, corresponding to OB and AGB stellar winds. The timestep-independent integration method matches the baseline to within $2\%$ (set by the accuracy of the implemented fits), whereas coarse discretization of the instantaneous rates suppresses mass loss, primarily due to unresolved early stellar evolution.
• Figure 3: Galaxy stellar mass function (GSMF) measurements from the calibrated CRK-HACC simulation at redshifts $z = 0$, $1$, and $2$, with red shaded regions indicating Poisson uncertainties. Gray shaded regions mark the resolution limit of the CRK-HACC results. Observational datasets bernardi2017driver2022thorne2021weaver2023 and empirical predictions from UniverseMachine (behroozi2019; the calibration target) are included for reference. Results from selected cosmological hydrodynamic simulations --- EAGLE furlong2015, SIMBAdave2019, COLIBREchaikin2025, Illustris-TNGpillepich2018b, FLAMINGOschaye2023flamingo, and MillenniumTNG pakmor2023 --- are also shown for qualitative comparison. The CRK-HACC GSMFs closely match both the observational constraints and the UniverseMachine predictions across all presented redshifts, including at high redshift where the spread among different simulation predictions becomes more pronounced.
• Figure 4: Median cluster gas density profiles from the fiducial CRK-HACC simulation at $z\!= \!0$ for halos with $M_{500c} > 3\times10^{14}\,M_\odot$. The shaded region shows the 16th-84th percentile range. Observational measurements are shown from mcdonald2017 (calibration target), ghirardini2021, and lyskova2023. Results from recent cosmological simulations include TNG-Cluster lehle2024, MillenniumTNG pakmor2023, FLAMINGO braspenning2024, and the SIMBA and GADGET-X clusters from The Three Hundred project li2023. The calibrated CRK-HACC profile lies near the midpoint of both the observational and simulation ranges, indicating consistency with current constraints on intracluster gas structure.
• Figure 5: The large central panel shows a slice of the gas density field across the full simulation domain ($256\,h^{-1}\mathrm{Mpc}$ on a side) with a slice depth of $4\,h^{-1}\mathrm{Mpc}$. The left half of this panel is color-mapped by temperature (cold in blue, hot in red), while the right half is color-mapped by metallicity (metal-poor in blue, metal-rich in green). The panels along the left column present successive zooms into the most massive halo in the volume ($M_{200c} = 1.2\times 10^{15}\,h^{-1}M_\odot$), and the panels along the bottom row show progressive zooms into a second massive cluster ($M_{200c} = 1.4\times 10^{14}\,h^{-1}M_\odot$). The lower-left panel shows the stellar density field in the central region of the largest halo.