GALFORM: Semi-Analytic Galaxy Formation Model

• GALFORM is a semi-analytic model of galaxy formation and evolution that simulates the build-up of galaxies within hierarchically growing dark-matter haloes in a ΛCDM universe.
• It employs coupled prescriptions for gas cooling, disc and burst star formation, feedback, and chemical enrichment to evolve baryonic reservoirs drawn from N-body or Monte Carlo merger trees.
• Widely calibrated against multi-wavelength observations and simulations, GALFORM serves as both a predictive tool for galaxy statistics and a computational platform for survey-mock predictions.

GALFORM is a semi-analytic model of galaxy formation and evolution developed in the Durham framework to follow the build-up of galaxies within hierarchically growing dark-matter haloes in a $\Lambda$CDM universe (Vernon et al., 2014). It is implemented on dark-matter merger trees extracted either from large $N$-body simulations or from Monte Carlo trees, and evolves baryonic reservoirs through coupled prescriptions for gas cooling, disc formation, star formation, feedback, chemical enrichment, mergers, disc instabilities, stellar population synthesis, and dust radiative transfer (Lacey et al., 2015). Across its major releases, GALFORM has been used both as a physical forward model for galaxy populations and as a computational platform for clustering, line-emission, cold-gas, radio, and survey-mock predictions, while also serving as a testbed for emulator-based calibration and for direct comparison with hydrodynamical simulations such as EAGLE (Mitchell et al., 2017).

1. Framework and model architecture

GALFORM follows galaxy formation by grafting analytic baryonic prescriptions onto merger trees of dark-matter haloes (Gonzalez-Perez et al., 2013). In the WMAP7 implementation discussed by Gonzalez-Perez et al., the underlying MS-W7 $N$-body simulation has a box size of $500\,h^{-1}\,\mathrm{Mpc}$, $2160^3$ particles, and 61 snapshots from $z=127$ to $z=0$, with haloes identified using Friends-of-Friends and SUBFIND and merger trees built with a core-tracking algorithm (Gonzalez-Perez et al., 2013). In the Planck-Millennium implementation, the simulation volume is $(800\,\mathrm{Mpc})^3$ with $5040^3$ particles and 269 snapshots, and merger trees are built with the DHALOS algorithm (Huško et al., 2022).

A central structural feature of GALFORM is its reservoir-based treatment of baryons. Mitchell et al. write the model in terms of four reservoirs—diffuse halo gas, ISM gas, stars, and an ejected reservoir—evolved through global, linearised rate equations (Mitchell et al., 2017). In that formulation, cooling or infall transfers gas from the diffuse halo into the ISM on a timescale $\tau_{\rm infall}$, star formation consumes ISM on $N$0, supernova-driven outflows populate the ejected reservoir with mass-loading factor $N$1, and reincorporation returns ejected gas to the diffuse halo on $N$2 (Mitchell et al., 2017). This closed system enforces mass conservation and, in its full form, also evolves the total angular momentum and metal mass in each reservoir (Mitchell et al., 2017).

The model architecture has changed across releases. The version summarized in "Galaxy Formation Spanning Cosmic History" removes the earlier assumption of discrete halo “formation events,” so that halo virial quantities, cooling integrals and satellite orbits evolve smoothly at each timestep (Benson et al., 2010). Lacey et al. then describe a unified multi-wavelength version that combines a different initial mass function in quiescent star formation and in starbursts, AGN feedback suppressing gas cooling in massive haloes, an empirical star-formation law in galaxy discs based on molecular gas content, updated cosmology, a more accurate treatment of dynamical friction acting on satellite galaxies, and an updated stellar population model (Lacey et al., 2015). This suggests that GALFORM is best regarded not as a single frozen prescription but as a family of closely related semi-analytic realizations calibrated against different observational targets.

2. Core baryonic physics

Gas cooling in GALFORM is treated by assuming that baryons in collapsed haloes are shock-heated to the virial temperature and then cool radiatively onto the central galaxy (Lacey et al., 2015). In one widely used formulation, the local cooling time is

$N$3

with cooling and free-fall radii defining an accretion radius $N$4 and a cooling rate $N$5 (Lacey et al., 2015). Earlier and later formulations differ in detail: Bower et al. emphasize continuous time-available cooling rather than resets at halo-doubling events (Benson et al., 2010), whereas Mitchell et al. identify the GALFORM cooling module as using a notional hot-gas profile and computing a cooling radius by equating the local cooling time to the time since the halo last doubled its mass (Mitchell et al., 2017).

Star formation is split into quiescent disc mode and burst mode. In the molecular-gas-based disc law introduced in the Lagos et al. line of GALFORM development, the star-formation surface density is

$N$6

with molecular fraction tied to mid-plane pressure through the Blitz and Rosolowsky relation (Lagos et al., 2012). In the Lacey et al. unified model, the global disc law is written as

$N$7

with $N$8 in that calibration (Lacey et al., 2015). By contrast, older implementations used $N$9 or $N$0 forms [(Vernon et al., 2014); (Benson et al., 2010)]. Starbursts are triggered by mergers or disc instabilities, and typically consume gas on a dynamical-time-related timescale such as

$N$1

so that $N$2 (Lacey et al., 2015, Madar et al., 2024).

Supernova feedback is represented through a mass-loading law in circular velocity. A common GALFORM form is

$N$3

with the normalization and slope treated as free parameters in calibration studies (Elliott et al., 2021, Madar et al., 2024). In the Lacey et al. model, the standard values are $N$4 and $N$5, with ejected gas returning on $N$6 and $N$7 (Lacey et al., 2015). The 2016 GSMF study further distinguishes $N$8 and $N$9, showing that the GSMF evolution strongly constrains quiescent-disc feedback while leaving burst loading more weakly restricted (Rodrigues et al., 2016).

AGN feedback in GALFORM is generally implemented as radio-mode suppression of cooling in massive haloes (Lacey et al., 2015). In one form, cooling is switched off when

$500\,h^{-1}\,\mathrm{Mpc}$0

with black-hole growth occurring during starbursts and in hot haloes (Lacey et al., 2015). Mitchell et al. stress that this is a bimodal “all or nothing” suppression: once the quasi-hydrostatic criterion is met and AGN heating exceeds the cooling luminosity, infall from the diffuse halo is completely shut down, but no hot gas is ejected from the halo (Mitchell et al., 2017).

Chemical enrichment and recycling are treated self-consistently. The model tracks recycled fraction $500\,h^{-1}\,\mathrm{Mpc}$1 and yield $500\,h^{-1}\,\mathrm{Mpc}$2, which depend on the IMF, and follows the exchange of metals between cold gas, hot gas, stars, and ejected reservoirs (Lacey et al., 2015). Some versions adopt instantaneous recycling (Lacey et al., 2015), while the 2010 model includes a non-instantaneous treatment following stellar lifetimes and yields (Benson et al., 2010).

3. Initial mass functions, dynamics, and satellite evolution

A defining aspect of several GALFORM variants is the use of different initial mass functions in quiescent and burst star formation. In the IMF definition

$500\,h^{-1}\,\mathrm{Mpc}$3

the Gon14 model adopts a universal Kennicutt IMF with slope $500\,h^{-1}\,\mathrm{Mpc}$4 for $500\,h^{-1}\,\mathrm{Mpc}$5 and $500\,h^{-1}\,\mathrm{Mpc}$6 for $500\,h^{-1}\,\mathrm{Mpc}$7, while Lac14 uses the same Kennicutt IMF in quiescent star formation but a single power law with top-heavy slope $500\,h^{-1}\,\mathrm{Mpc}$8 in starbursts (Campbell et al., 2014). Lacey et al. similarly employ a Kennicutt IMF in discs and a mildly top-heavy burst IMF with $500\,h^{-1}\,\mathrm{Mpc}$9, concluding that the observed number counts and redshift distribution of sub-millimetre galaxies require a burst IMF that is somewhat top-heavy (Lacey et al., 2015).

Galaxy mergers and disc instabilities provide the main channels for spheroid formation and burst triggering. Major mergers destroy discs and trigger starbursts when the baryonic mass ratio exceeds $2160^3$0, while disc instabilities are triggered when

$2160^3$1

or in the closely related Efstathiou-form expression

$2160^3$2

(Lacey et al., 2015, Huško et al., 2022, Madar et al., 2024). Huško et al. show that in the Planck-Millennium GALFORM run, mergers dominate in building up high-mass and low-mass spheroids, while disc instabilities and their associated starbursts dominate for intermediate-mass spheroids at $2160^3$3 (Huško et al., 2022).

The treatment of satellites is one of the most consequential dynamical components. In standard GALFORM, a galaxy entering another FoF halo becomes a satellite and remains so thereafter, with dynamical-friction merger times assigned analytically (Lacerna et al., 2017). In the default satellite-merger treatment discussed by Campbell et al., an analytic dynamical friction timescale $2160^3$4 is immediately assigned when a galaxy becomes a satellite, without tracking surviving subhaloes (Campbell et al., 2014). Their subhalo dynamical friction (“sdf”) scheme instead keeps satellites linked to resolved subhaloes until the subhalo drops below the simulation resolution; only then is $2160^3$5 computed from the last-known most-bound particle in an NFW host halo potential (Campbell et al., 2014). This delays mergers, yields a more extended radial distribution of satellites, and reduces clustering amplitude on scales $2160^3$6 by up to $2160^3$7–$2160^3$8 percent while leaving large-scale clustering almost unchanged (Campbell et al., 2014).

Environmental quenching prescriptions also differ between GALFORM variants. In the conformity study, GP14 removes all hot gas from a satellite instantaneously upon infall, whereas GRP uses gradual ram-pressure stripping, preserving some hot reservoir for longer and slowing quenching (Lacerna et al., 2017). This difference materially affects one-halo conformity signals, with GP14 over-quenching satellites and GRP yielding a stronger one-halo signal (Lacerna et al., 2017).

4. Angular momentum, sizes, SPS, dust, and line-emission modules

GALFORM assigns cooled gas the specific angular momentum of the dark matter halo and assumes that newly formed stars inherit the same specific angular momentum as the total gas disc:

$2160^3$9

Mitchell et al. show, by comparison to EAGLE, that this simplifying assumption overestimates the stellar specific angular momenta and hence galaxy sizes, because in EAGLE stars form preferentially out of low-$z=127$0 gas in the ISM (Mitchell et al., 2017). They identify galaxy sizes too large by up to $z=127$1 at $z=127$2 and propose replacing $z=127$3 with a star-forming-gas-weighted $z=127$4 (Mitchell et al., 2017). This criticism aligns with broader comparisons in which GALFORM overpredicts low-mass galaxy ISM fractions at $z=127$5 as a consequence of its large assumed stellar specific angular momentum and low gas surface density (Mitchell et al., 2017).

Stellar population synthesis and dust treatment are integral to GALFORM’s predictive scope. Different releases use different SPS libraries: Gon14 employs a private 1999 release of Bruzual and Charlot models, Lac14 uses Maraston (2005), and the unified multi-wavelength model uses Maraston (2005) including TP-AGB stars [(Campbell et al., 2014); (Lacey et al., 2015)]. The sensitivity study of Gonzalez-Perez et al. runs GALFORM with seven different SPS models and finds that rest-frame UV and optical luminosity functions are insensitive to SPS choice, whereas the rest-frame NIR luminosity function depends strongly on the treatment of the TP-AGB phase, with differences larger than a factor of 2 for galaxies brighter than $z=127$6 (Gonzalez-Perez et al., 2013). That result has direct methodological implications for stellar-mass estimation and for comparisons between model and data.

Dust attenuation and emission are treated with a two-phase radiative-transfer framework. Lacey et al. describe a two-phase ISM with molecular clouds and diffuse dust, age-dependent attenuation, and energy-balance dust emission with separate temperatures for clouds and diffuse material (Lacey et al., 2015). In clustering studies, GALFORM broad-band magnitudes are “observed” through the same filter sets used in surveys, with the model’s internal radiative-transfer dust attenuation included before SED fitting (Campbell et al., 2014). Campbell et al. then show that SED-fitted stellar masses can differ sharply from true model masses, with strong internal dust attenuation driving estimated masses well below the true value by up to an order of magnitude in dusty, high-mass systems (Campbell et al., 2014).

GALFORM has also been extended to predict nebular emission lines. In the star-forming-galaxy line model, GALFORM computes the hydrogen-ionizing photon rate

$z=127$7

and then assigns each emission line a luminosity

$z=127$8

where $z=127$9 is interpolated from a pre-computed Gutkin et al. CLOUDY grid as a function of ionisation parameter, metallicity, and hydrogen density (Baugh et al., 2021). Empirical relations from Kashino and Inoue then link $z=0$0 and $z=0$1 to galaxy stellar mass, gas-phase metallicity, and $z=0$2 (Baugh et al., 2021). This module reproduces the local BPT locus, predicts evolution in line-ratio diagrams, and indicates that high-redshift star-forming galaxies have gas densities and ionisation parameters approximately $z=0$3–$z=0$4 times higher than local systems in the relevant selections (Baugh et al., 2021).

Two further extensions illustrate the model’s modularity. GALFORM has been coupled to UCL_PDR to predict CO line emission from the rotational transitions $z=0$5–$z=0$6 to $z=0$7–$z=0$8, using predicted $z=0$9, gas metallicity, size, and radiation fields as inputs to a photon-dominated-region calculation (Lagos et al., 2012). It has also been combined with a Monte Carlo Ly$(800\,\mathrm{Mpc})^3$0 radiative-transfer code using Shell and Wind outflow geometries to predict Ly$(800\,\mathrm{Mpc})^3$1 escape fractions, line profiles, equivalent widths, and luminosity functions from $(800\,\mathrm{Mpc})^3$2 to $(800\,\mathrm{Mpc})^3$3 (Orsi et al., 2011).

5. Calibration, history matching, and deep-learning emulation

The computational cost of GALFORM has made statistical calibration a major research topic in its own right. Vernon, Goldstein, and Bower describe GALFORM as a code with many input parameters and simulations that take approximately $(800\,\mathrm{Mpc})^3$4 on one modern CPU, making brute-force comparison with observations infeasible (Vernon et al., 2014). Their solution is Bayesian emulation within an iterative history-matching strategy. For each output $(800\,\mathrm{Mpc})^3$5 they write an emulator of the form

$(800\,\mathrm{Mpc})^3$6

combine this with an observation model $(800\,\mathrm{Mpc})^3$7, and define an implausibility measure

$(800\,\mathrm{Mpc})^3$8

Across four main waves plus a final confirmation wave, the surviving region shrinks from the full 17-dimensional cube to $(800\,\mathrm{Mpc})^3$9 of the original volume (Vernon et al., 2014).

This framework was extended in the GSMF-evolution study to a 20-dimensional parameter space, where Bayesian emulation and history matching were used to constrain the model with stellar-mass functions from $5040^3$0 to $5040^3$1 (Rodrigues et al., 2016). That analysis finds that the GSMF strongly constrains parameters controlling quiescent star formation, stellar feedback, AGN feedback, and disc instability thresholds, but local data alone do not usually select models that match the GSMF evolution well (Rodrigues et al., 2016). Principal-component analysis then yields explicit relations among the most constrained parameter combinations, particularly those linking the supernova-feedback normalization in discs, the velocity dependence of mass loading, and the AGN threshold parameter (Rodrigues et al., 2016).

Deep-learning emulation is the next methodological step. In the 2021 calibration study, the updated GALFORM model is emulated with an ensemble of dense feedforward networks trained on $5040^3$2 GALFORM runs sampled via a Latin hypercube (Elliott et al., 2021). The architecture uses 2 hidden layers of 512 neurons with sigmoid activations, MAE loss, AMSGRAD followed by RMSprop fine-tuning, and an ensemble of 10 networks (Elliott et al., 2021). The resulting emulator achieves an out-of-sample mean absolute error at the knee of the K-band luminosity function of $5040^3$3 dex with less than $5040^3$4 model evaluations, and scaling approximately $5040^3$5 in training size (Elliott et al., 2021). Sensitivity analysis shows that more than $5040^3$6 of the variance in key observables arises from 2–3 parameters, with the faint luminosity-function end dominated by supernova-feedback parameters and the bright end primarily by $5040^3$7 (Elliott et al., 2021).

A directly related strategy underlies the 2024 H$5040^3$8-emitter forecasts. Madar et al. train a deep-learning emulator on 3000 GALFORM runs across an eleven-dimensional parameter space, using the local $5040^3$9- and $\tau_{\rm infall}$0-band luminosity functions together with H$\tau_{\rm infall}$1 redshift distributions as training observables (Madar et al., 2024). The network is a fully connected feed-forward ensemble of 5 identical networks with six hidden layers of 512 neurons and Leaky-ReLU activation, optimized with Adam+AMSGrad and then RMSprop (Madar et al., 2024). MCMC calibration with heuristic weights reveals tensions between local luminosity functions and higher-redshift H$\tau_{\rm infall}$2 counts, but finds improved fits to the H$\tau_{\rm infall}$3 number counts while maintaining appropriate predictions for the local universe luminosity function (Madar et al., 2024).

6. Empirical tests, survey applications, and known limitations

GALFORM has been tested against a wide range of observables. The unified multi-wavelength model is reported to reproduce the observed evolution of the $\tau_{\rm infall}$4-band luminosity function and stellar-mass function, the number counts and redshift distribution of $\tau_{\rm infall}$5-selected sub-mm galaxies, the $\tau_{\rm infall}$6 optical and near-IR luminosity functions, the HI mass function, the fraction of early-type galaxies, the Tully–Fisher relation, metallicity-luminosity and size-luminosity relations at $\tau_{\rm infall}$7, far-IR number counts, and far-UV luminosity functions at $\tau_{\rm infall}$8–$\tau_{\rm infall}$9 (Lacey et al., 2015). At the same time, Lacey et al. note residual discrepancies: sizes of low-luminosity discs and spheroids are too large, stellar metallicities of low-mass early types are underpredicted, colour bimodality is too strong, the low-mass end of the high-$N$00 stellar-mass function is too steep, and the cosmic star-formation rate density is somewhat low at $N$01 (Lacey et al., 2015).

Clustering tests have become a particularly stringent application. Campbell et al. use GALFORM to predict the two-point correlation function of galaxy clustering as a function of stellar mass, embedding galaxies in the MS-W7 simulation and measuring $N$02 in 30 logarithmic bins over $N$03 (Campbell et al., 2014). They define the projected correlation function

$N$04

and show that comparing theory to observation requires applying to the model the same SED fitting and abundance-matching steps used in observational analyses (Campbell et al., 2014). With the sdf merger scheme, SED fitting, and abundance matching, agreement improves for SDSS and GAMA on small scales, but discrepancies remain for some stellar-mass ranges and separation scales, and VIPERS still exposes model-dependent failures (Campbell et al., 2014).

GALFORM has also been used to study galactic conformity. In GP14 and GRP, a primary galaxy is selected through an isolation criterion, quenched fractions are measured as

$N$05

and conformity is defined as

$N$06

The two-halo conformity signal in GALFORM is only at the few-percent level beyond $N$07, and largely disappears when satellites are rigorously purged from the primary sample (Lacerna et al., 2017). This weak long-range signal contrasts with earlier strong claims from SDSS and points to contamination by misclassified satellites and halo-mass scatter at fixed stellar mass as major confounding factors (Lacerna et al., 2017).

Survey-mock construction is another mature application. Izquierdo-Villalba et al. build a “digital twin” of a luminous red galaxy spectroscopic survey by running GALFORM on the high-resolution Planck-Millennium simulation and grafting its halo occupation predictions onto one thousand larger GLAM simulations (Hernández-Aguayo et al., 2020). The central and satellite occupations $N$08 and $N$09 are tabulated from GALFORM and used to populate GLAM haloes, with centrals at halo centres and satellites distributed according to an NFW profile (Hernández-Aguayo et al., 2020). This produces a non-standard HOD for DESI-like LRGs: $N$10 never reaches unity, peaks at $N$11 at $N$12, and declines at higher masses, while the satellite fraction drops from approximately $N$13 to $N$14 over $N$15 (Hernández-Aguayo et al., 2020). The resulting mocks reproduce GALFORM clustering to few-percent accuracy, match recent photometric clustering measurements of DESI-like LRGs, and yield covariance matrices plus forecasts for the linear-growth rate and BAO distances (Hernández-Aguayo et al., 2020).

The H$N$16-emitter forecasts provide an analogous application for future redshift surveys. With the recalibrated model, a Euclid-like survey with flux limit $N$17 over $N$18 yields predicted surface densities of $N$19–$N$20 in the 10–90 percentile range in the detailed report, while the abstract reports $N$21–$N$22; for a Roman-like survey with flux limit $N$23 and $N$24, the corresponding ranges are $N$25–$N$26 in the detailed report and $N$27–$N$28 in the abstract (Madar et al., 2024). This suggests that the precise forecast interval depends on which summary statistic is quoted, but the common conclusion is that GALFORM can be calibrated directly to H$N$29 number-count data and then used to predict effective bias and source densities for planned surveys.

Several limitations recur across independent studies. Comparison with EAGLE identifies overly rapid baryon cycling, an artificial dependence of infall and AGN feedback on halo mass-doubling events, and an angular-momentum prescription that overestimates galaxy sizes (Mitchell et al., 2017). The CSED analysis finds that GALFORM predicts the CSED for $N$30 in good agreement with observations but becomes increasingly poor toward $N$31, where the model CSED is $N$32 per cent fainter; the integrated galactic light is then $N$33 per cent fainter, consistent with the model underpredicting the cosmic star formation history (Andrews et al., 2017). The spheroid-assembly study finds ex-situ fractions at $N$34 that are $N$35–$N$36 higher than in EAGLE, TNG, and EMERGE, a high-mass ex-situ relation almost independent of redshift, and a peak pseudobulge fraction of $N$37 at $N$38, likely too high relative to observations (Huško et al., 2022). In radio applications, GALFORM coupled to a dynamo model reproduces the observed $N$39 radio luminosity function up to $N$40 only when turbulent speed is correlated with star-formation rate; models with constant or purely redshift-driven turbulence perform poorly, and the formalism appears inadequate once high-redshift galaxy structure and turbulence depart strongly from their low-$N$41 analogues (Jose et al., 2024).

A plausible implication is that GALFORM’s enduring strength lies in its modularity and calibration efficiency rather than in the literal finality of any single prescription. The model has repeatedly been extended—through new star-formation laws, revised satellite dynamics, top-heavy burst IMFs, SPS updates, emission-line grids, PDR coupling, Ly$N$42 radiative transfer, Bayesian emulators, and deep-learning surrogates—while preserving a common reservoir-based semi-analytic backbone (Lacey et al., 2015). For that reason, GALFORM remains both a scientific model of galaxy formation physics and an evolving inference framework for confronting multi-wavelength, multi-redshift survey data with physically structured but computationally tractable predictions.