COLIBRE: Cold ISM & Galaxy Formation Simulations

• COLIBRE simulations are a suite of cosmological hydrodynamical models that explicitly resolve a cold, multiphase ISM with live dust evolution and non-equilibrium chemistry.
• They employ innovative methodologies such as dark-matter supersampling and forward modelling to reproduce observed properties like stellar mass functions, galaxy sizes, and black-hole masses.
• The integrated baryonic physics, including multi-channel stellar feedback and AGN implementations, enables studies of galaxy evolution, star formation laws, and quenching across cosmic time.

COLIBRE is a suite of cosmological hydrodynamical simulations and an associated galaxy-formation model developed to follow galaxy assembly from very high redshift to the present day while explicitly modelling a cold, multiphase interstellar medium (ISM), live dust evolution, non-equilibrium chemistry, and calibrated stellar and black-hole feedback within representative cosmological volumes. The acronym denotes COld ISM and Better REsolution. A defining design choice is the removal of an imposed ISM pressure floor, allowing gas to cool to $\sim 10$ K; another is dark-matter supersampling by a factor of 4, which reduces spurious dark-matter-to-stellar heating by making dark-matter and baryonic particle masses comparable. The suite is implemented in SWIFT with the SPHENIX hydrodynamics scheme and was calibrated primarily to reproduce the $z\approx 0$ stellar mass function, galaxy sizes, and black-hole masses, after which it has been used to study stellar masses, star formation, metallicities, dust, morphologies, quenching, tidal stripping, and luminosity functions from the far-UV to the submillimetre (Schaye et al., 28 Aug 2025).

1. Numerical design and simulation architecture

COLIBRE adopts a “wedding cake” strategy in box size and mass resolution. The suite spans baryonic particle masses of $\sim 10^5$, $10^6$, and $10^7\,M_\odot$, with maximum box sizes of 50, 200, and 400 cMpc, respectively. The flagship large-volume runs include L200m6 and L400m7, and the two largest runs contain $5\times 3008^3$ particles, i.e. about $1.4\times 10^{11}$ particles (Schaye et al., 28 Aug 2025).

Resolution	Baryonic particle mass	Largest box
m5	$2.30\times 10^5\,M_\odot$	50 cMpc
m6	$1.84\times 10^6\,M_\odot$	200 cMpc
m7	$1.47\times 10^7\,M_\odot$	400 cMpc

Initial conditions are generated with monofonIC/MONOFONIC at $z\approx 0$0, using second-order Lagrangian perturbation theory and partially fixed amplitudes to reduce cosmic variance. The adopted cosmology is $z\approx 0$1, $z\approx 0$2, $z\approx 0$3, $z\approx 0$4, $z\approx 0$5, with one massive neutrino of $z\approx 0$6 eV treated semi-linearly rather than with neutrino particles (Schaye et al., 28 Aug 2025).

The numerical architecture combines a particle-mesh solver for long-range gravity with a fourth-order fast multipole method for short-range gravity, and uses the SPHENIX density-energy SPH formulation for gas dynamics. A central numerical innovation is the 4:1 dark-matter-to-baryon particle-number ratio. Because this brings dark-matter and baryonic particle masses close to parity, it suppresses spurious two-body energy transfer from dark matter to stars; the morphology and size studies identify this as important for retaining ordered rotation and avoiding artificially inflated stellar systems (Moreno et al., 3 Apr 2026).

2. Baryonic physics and subgrid model

The COLIBRE baryonic model combines radiative cooling and heating, non-equilibrium H/He chemistry, dust formation and evolution, star formation, stellar mass loss and enrichment, turbulent diffusion, early stellar feedback, supernova feedback, black-hole growth, and AGN feedback within a unified framework (Schaye et al., 28 Aug 2025).

Cooling is computed with hybrid-chimes/CHIMES, down to $z\approx 0$7 K. Hydrogen and helium are evolved out of equilibrium, their free-electron contribution enters metal-line cooling, and the network is coupled to shielding, a metagalactic UV/X-ray background, a local interstellar radiation field, and cosmic rays. The chemical and thermal model is dust-coupled rather than appended in post-processing: COLIBRE tracks three grain species and two grain sizes, with growth, destruction, shattering, coagulation, sputtering, astration, turbulent diffusion, and feedback-driven destruction all represented on the fly (Schaye et al., 28 Aug 2025).

Star formation follows a Schmidt-law-like prescription with a fixed efficiency per free-fall time,

$z\approx 0$8

applied only to gas that satisfies a gravitational instability criterion at the resolution scale. In this sense, the simulation imposes neither a molecular Kennicutt–Schmidt relation nor an effective ISM equation of state; instead, the resolved H I and H$z\approx 0$9 star-formation laws emerge from the interaction of cooling, chemistry, dust, and feedback (Lagos et al., 12 Dec 2025).

Stellar feedback is explicitly multi-channel. Before the first core-collapse supernovae, COLIBRE includes non-explosive pre-supernova (NEPS) feedback from stellar winds, radiation pressure, and $\sim 10^5$0-region photoheating. The NEPS implementation derives age- and metallicity-dependent energy and momentum budgets from BPASS and, in isolated-disk tests, regulates star formation chiefly through $\sim 10^5$1-region heating, improves convergence, and pre-processes the ISM for later supernova explosions (Benítez-Llambay et al., 29 Sep 2025). Core-collapse supernova feedback itself is stochastic thermal+kinetic feedback with a fiducial kinetic fraction $\sim 10^5$2, low-velocity $\sim 10^5$3 kicks, a density-dependent thermal heating temperature, and a pressure-dependent energy fraction $\sim 10^5$4; the calibration study found that this environmental dependence was required to reproduce both galaxy masses and galaxy sizes simultaneously (Chaikin et al., 4 Sep 2025).

Black holes are seeded in sufficiently massive haloes and grow through modified Bondi–Hoyle–Lyttleton accretion with turbulence and vorticity corrections, capped at 100 times the Eddington rate. COLIBRE provides two AGN implementations: a fiducial thermal model in which stored energy is injected thermally with a black-hole-mass-dependent heating temperature, and a hybrid model that adds spin evolution, accretion-disc state changes, thermal winds, and kinetic jets. The hybrid model is intended to probe uncertainties in AGN coupling, not merely to provide a numerically different variant (Schaye et al., 28 Aug 2025).

3. Calibration strategy and analysis ecosystem

COLIBRE is not presented as an ab initio calculation without calibration. The subgrid model was explicitly tuned to match the observed $\sim 10^5$5 galaxy stellar mass function, the size–stellar mass relation, and the black-hole mass–stellar mass relation in massive galaxies. The principal calibration campaign used Latin hypercubes of approximately 200 cosmological m7 runs in $\sim 10^5$6 volumes, with Gaussian-process emulators trained on the simulation outputs and then sampled with MCMC. Simpler feedback models could reproduce either the GSMF or the size relation in isolation, but the full COLIBRE model with density-dependent thermal supernova heating and pressure-dependent supernova energy achieved the best simultaneous fit, with $\sim 10^5$7 in the calibration study (Chaikin et al., 4 Sep 2025).

Galaxy and subhalo identification in COLIBRE commonly uses HBT-HERONS, a history-based subhalo finder and merger-tree builder that is particularly robust for heavily stripped satellites. Several analyses measure galaxy properties with SOAP and adopt fiducial physical apertures such as 50 pkpc for stellar masses and star-formation rates. This analysis infrastructure matters because a number of COLIBRE results depend on distinguishing central from satellite evolution, following progenitor histories, or measuring apertures in ways that emulate observations (Chandro-Gómez et al., 29 May 2026).

A distinctive feature of the project is extensive forward modelling into observable space. High-redshift UV luminosity functions are obtained by post-processing COLIBRE galaxies with SKIRT, using dust distributions predicted directly by the simulation, BPASS templates for evolved stars, and the TODDLERS library for star-forming regions younger than 10 Myr. For the UVLF work, SKIRT is run in ExtinctionOnly mode, with UV magnitudes measured at 1500 Å in a 50 kpc proper aperture (Lu et al., 7 May 2026). At $\sim 10^5$8, a related COLIBRE+SKIRT pipeline predicts luminosity functions from the far-UV to $\sim 10^5$9, using dust attenuation and thermal re-emission from the live simulation dust field rather than an assumed fixed dust-to-metal ratio (Lu et al., 3 May 2026).

Observable definitions are correspondingly explicit. For example, the gas-phase mass–metallicity study measures oxygen abundance from cool, dense star-forming gas with $10^6$0 and $10^6$1, uses a mass-weighted oxygen-to-hydrogen ratio, and measures metallicity in a 3 kpc aperture while stellar masses and SFRs are measured within 50 kpc to mimic common observational practice (Sharda et al., 24 Jun 2026).

4. Predictions for galaxy populations across cosmic time

A central COLIBRE result is that the galaxy stellar mass function evolves smoothly from $10^6$2 to $10^6$3, with broad agreement with observations over $10^6$4 and maximum systematic deviations of about $10^6$5–$10^6$6 dex at $10^6$7. The same study finds good agreement for the star-forming main sequence, cosmic star-formation-rate density, stellar mass density, and the quenched fraction, while arguing that neither a redshift-dependent star-formation efficiency, nor a variable IMF, nor non-$10^6$8CDM cosmology is required to reproduce JWST-era stellar masses and SFRs (Chaikin et al., 9 Sep 2025).

The gas-phase mass–metallicity relation is another major COLIBRE benchmark. In the simulations, the MZR is already established by $10^6$9, shows little evolution until $10^7\,M_\odot$0, and becomes shallower toward low redshift. Its low-mass slope is controlled mainly by core-collapse supernova feedback, while the high-mass turnover and downturn are largely governed by AGN feedback. Variations in star-formation efficiency or oxygen depletion onto dust grains affect the MZR more weakly, although ignoring oxygen depletion can shift the high-mass relation upward by $10^7\,M_\odot$1 dex and in some cases by $10^7\,M_\odot$2 dex (Sharda et al., 24 Jun 2026).

On kiloparsec scales, COLIBRE predicts resolved H I, H$10^7\,M_\odot$3, and total-neutral-gas Kennicutt–Schmidt relations without explicitly imposing them. At $10^7\,M_\odot$4, the resolved relations and their scatter agree with observations, while from $10^7\,M_\odot$5 to $10^7\,M_\odot$6 the molecular depletion time decreases by a factor of $10^7\,M_\odot$7, primarily because lower gas-phase metallicities at high redshift shift the cold ISM toward a more atomic state at fixed star-formation rate (Lagos et al., 12 Dec 2025).

High-redshift UV luminosity functions show a more qualified outcome. COLIBRE reproduces the observed evolution of the stellar mass function up to $10^7\,M_\odot$8, but its dust-attenuated UVLFs are too faint at the bright end: at a number density of $10^7\,M_\odot$9, the brightest galaxies are underluminous by $5\times 3008^3$0 mag at $5\times 3008^3$1 and $5\times 3008^3$2 mag at $5\times 3008^3$3. Removing dust improves agreement at $5\times 3008^3$4, but even dust-free galaxies remain too faint at $5\times 3008^3$5, leading that study to suggest that additional physics may be required for the earliest UV-bright systems, with a top-heavy IMF offered as a leading possibility (Lu et al., 7 May 2026).

5. Structure, morphology, satellites, and quenching

COLIBRE has been used extensively to study structural scaling relations. At $5\times 3008^3$6, it reproduces observed size–mass relations over $5\times 3008^3$7 for multiple size definitions, and recovers the morphology-dependent segregation of galaxies in both the size–mass and stellar-specific-angular-momentum planes. Star-forming discs and quenched spheroids occupy distinct, approximately parallel $5\times 3008^3$8–$5\times 3008^3$9 sequences, with an offset of roughly $1.4\times 10^{11}$0–$1.4\times 10^{11}$1 dex, and agreement with observed angular momentum content remains good out to $1.4\times 10^{11}$2 (Ludlow et al., 27 Mar 2026).

The present-day morphology study argues that COLIBRE’s cold-ISM treatment and dark-matter supersampling directly address two numerical problems that historically affected representative-volume simulations: over-pressurised gas discs and spurious stellar heating by massive dark-matter particles. Morphology was not a calibration target, yet median morphology–mass relations are well converged across m5, m6, and m7. The simulations predict that galaxies with stellar masses of $1.4\times 10^{11}$3 are the most rotationally dominated, that morphology correlates only weakly with halo properties at fixed stellar mass, and that stronger trends are instead seen with internal gas content, SFR, stellar age, dust, black-hole mass, ex-situ fraction, and metallicity (Moreno et al., 3 Apr 2026).

High-redshift massive quenched galaxies constitute another major application. In COLIBRE, MQGs are selected with $1.4\times 10^{11}$4 and $1.4\times 10^{11}$5. The simulations find number densities and stellar mass functions in broad agreement with current observations once observational uncertainties are forward-modelled, and identify AGN feedback as the primary quenching mechanism. MQGs host more massive black holes and have higher star-formation efficiencies than equally massive star-forming systems; the origin of these differences is traced to overdense environments on $1.4\times 10^{11}$6 and $1.4\times 10^{11}$7 scales before quenching, where gas inflows, black-hole accretion, and feedback power are enhanced (Chandro-Gómez et al., 18 Dec 2025). A follow-up comparison between the thermal and hybrid AGN models shows that thermal AGN feedback quenches more efficiently at $1.4\times 10^{11}$8, whereas the hybrid model is delayed by slower early black-hole growth and the longer timescale over which jets affect host-galaxy gas (Chandro-Gómez et al., 29 May 2026).

Satellite evolution in COLIBRE is tracked with unusual fidelity by HBT-HERONS. The satellite stripping study identifies a universal tidal track between the remaining stellar mass fraction and the remaining subhalo mass fraction, with a characteristic transition at $1.4\times 10^{11}$9 and median late-stage slope $2.30\times 10^5\,M_\odot$0. Most of the early stripping is dark matter, while substantial stellar stripping begins only after the bound subhalo mass fraction falls below that threshold. The same work uses COLIBRE to explain the paucity of long-lived “orphan” galaxies in hydrodynamical simulations and predicts a substantial population of dark-matter-deficient galaxies, with abundance peaking at $2.30\times 10^5\,M_\odot$1 (He et al., 3 Apr 2026).

6. Limitations, tensions, and interpretive caveats

COLIBRE’s scope is broad, but its limitations are explicitly stated in the project papers. Molecular clouds are not resolved; several processes remain subgrid and calibrated rather than first-principles predictions; black-hole dynamics rely on approximate repositioning; the hydrodynamical runs do not include full radiative transfer; magnetic fields are absent; and some observational comparisons remain primarily in theory space unless forward-modelled. Numerical convergence is strong for many observables, but weaker for molecular gas fractions at low mass, dust masses and grain sizes in some regimes, and AGN-sensitive quantities because central gas densities increase with resolution (Schaye et al., 28 Aug 2025).

Several specific tensions have emerged. In local luminosity functions, COLIBRE+SKIRT matches data from the far-UV to $2.30\times 10^5\,M_\odot$2 and from 70 to $2.30\times 10^5\,M_\odot$3, but underpredicts the bright end of the mid-infrared luminosity functions, with the discrepancy growing from about $2.30\times 10^5\,M_\odot$4–$2.30\times 10^5\,M_\odot$5 dex at $2.30\times 10^5\,M_\odot$6 to $2.30\times 10^5\,M_\odot$7 dex at $2.30\times 10^5\,M_\odot$8 at a number density of $2.30\times 10^5\,M_\odot$9 (Lu et al., 3 May 2026). At high redshift, the brightest UV-selected galaxies remain too faint, especially by $1.84\times 10^6\,M_\odot$0, even after dust attenuation is removed and observational uncertainties are considered (Lu et al., 7 May 2026). Structural comparisons also show that simulated galaxies with $1.84\times 10^6\,M_\odot$1 are somewhat smaller than observed at $1.84\times 10^6\,M_\odot$2; the size study argues that neglected dust attenuation could increase apparent sizes by $1.84\times 10^6\,M_\odot$3–$1.84\times 10^6\,M_\odot$4 dex at $1.84\times 10^6\,M_\odot$5, especially where the discrepancy is strongest (Ludlow et al., 27 Mar 2026).

These caveats define COLIBRE’s current scientific status. It is neither a purely phenomenological post-processing framework nor a parameter-free theory. Rather, it is a calibrated cosmological simulation programme in which a cold multiphase ISM, dust-coupled thermochemistry, dark-matter supersampling, and forward-modelled observables are combined to make numerically converged and observationally testable predictions across many sectors of galaxy formation. The suite’s main significance lies in showing that many relations traditionally treated separately—stellar mass growth, gas metallicity, resolved star formation laws, dust-obscured emission, angular momentum, morphology, quenching, and tidal stripping—can be analysed within a single hydrodynamical framework whose strengths and residual tensions are both quantitatively exposed (Schaye et al., 28 Aug 2025).