FirstLight: Simulating Cosmic Dawn Galaxies

• FirstLight is a suite of hydrodynamical zoom-in simulations that model primeval galaxies during cosmic dawn and reionization using adaptive mesh refinement and detailed subgrid physics.
• It predicts key observables including UV luminosities, bursty star formation histories, nebular emissions, and dust-affected spectral energy distributions, linking simulations to JWST-era data.
• Key findings reveal redshift-dependent star formation efficiency, evolving mass-metallicity relations, and distinctive galaxy compactness that elucidate the growth and structure of early galaxies.

FirstLight is a suite of cosmological zoom-in hydrodynamical simulations constructed to follow the formation and evolution of primeval galaxies during cosmic dawn and the epoch of reionization, typically from $z \gtrsim 15$ down to $z \simeq 5$–6. Across the FirstLight series and its later extensions, the project couples high-resolution adaptive-mesh simulations to post-processing for stellar populations, nebular emission, dust radiative transfer, and synthetic observables, with the explicit aim of building statistically meaningful, mass-limited samples of early galaxies and predicting their UV luminosities, star-formation histories, metallicities, line emission, dust attenuation, morphologies, and sizes (Ceverino et al., 2017, Ceverino et al., 2021).

1. Project scope and development

FirstLight was introduced as a large database of high-resolution zoom-in simulations of galaxy formation around the epoch of reionization. The initial project design selected halos from parent boxes of 10, 20, and 40 $h^{-1}\,\mathrm{Mpc}$, with a main halo catalog of 978 objects across all boxes and cosmologies; the first published tests analyzed 15 zoom-in halos with $M_h \approx 10^9$–$10^{11}\,M_\odot$ at $z=5$ (Ceverino et al., 2017). Subsequent analyses focused on a mass-selected sample of 290 halos at $z=5$, yielding roughly 300 distinct galaxy histories over $z=5$–15, and later extensions expanded the database to 377 and 430 zoom-in simulations for studies of efficiency and synthetic size measurements (Ceverino et al., 2018, Ceverino et al., 2024, Ceverino et al., 5 Mar 2026).

The simulated galaxy population spans broad ranges in intrinsic properties. The database described in the early series covers rest-frame UV magnitudes $M_{\rm UV} \simeq -12$ to $-22$ and stellar masses $z \simeq 5$0–$z \simeq 5$1, while specific studies of burst histories quote $z \simeq 5$2–$z \simeq 5$3 (Ceverino et al., 2021, Ceverino et al., 2018). Later analyses of bright and massive systems extend to $z \simeq 5$4–$z \simeq 5$5 in selected subsamples (Nakazato et al., 2024, Ceverino et al., 5 Mar 2026).

The internal structure of the series is cumulative. FirstLight I established the UV luminosity function, stellar mass function, and basic halo–galaxy scaling relations at $z \simeq 5$6; FirstLight II quantified the highly bursty star-formation histories from $z \simeq 5$7 to 15; FirstLight III modeled rest-frame UV–optical spectral energy distributions; FirstLight IV added nebular-line physics and spatially resolved emission in sub-$z \simeq 5$8 systems at $z \simeq 5$9 (Ceverino et al., 2017, Ceverino et al., 2018, Ceverino et al., 2018, Ceverino et al., 2021). Later papers extended the framework to the mass–metallicity relation, [O III]-bright emitters, dust attenuation and infrared emission, clump formation, galaxy compactness, and accelerated size evolution (Langan et al., 2019, Nakazato et al., 2023, Mushtaq et al., 2023, Nakazato et al., 2024, Cataldi et al., 6 Oct 2025, Ceverino et al., 5 Mar 2026).

Series component	Primary focus
FirstLight I	UV luminosity function, SMF, scaling relations
FirstLight II	Bursty SFR histories and SF main sequence
FirstLight III	Rest-frame UV–optical SEDs
FirstLight IV	[O III], H$h^{-1}\,\mathrm{Mpc}$0, and spatially resolved nebular structure

This progression made FirstLight not merely a galaxy catalog, but an interpretive framework for JWST-era studies of reionization galaxies.

2. Numerical framework and subgrid physics

Most FirstLight analyses use the ART code, an adaptive mesh refinement Eulerian gravity-plus-hydrodynamics solver. In the standard database, the simulations adopt WMAP5 cosmology with $h^{-1}\,\mathrm{Mpc}$1, $h^{-1}\,\mathrm{Mpc}$2, $h^{-1}\,\mathrm{Mpc}$3, and $h^{-1}\,\mathrm{Mpc}$4; some later 40 $h^{-1}\,\mathrm{Mpc}$5 extensions use Planck13 or Planck 2014 cosmologies (Ceverino et al., 2018, Langan et al., 2019, Nakazato et al., 2024, Nakazato et al., 7 Feb 2026). Target halos are selected by $h^{-1}\,\mathrm{Mpc}$6 at $h^{-1}\,\mathrm{Mpc}$7, with thresholds that depend on box size, and are re-simulated at high resolution.

The canonical FirstLight resolution is $h^{-1}\,\mathrm{Mpc}$8 for dark matter, minimum star-particle mass $h^{-1}\,\mathrm{Mpc}$9, and maximum spatial resolution 8.7–17 proper pc; the more massive 40 $M_h \approx 10^9$0 samples typically use $M_h \approx 10^9$1, minimum stellar particle mass $M_h \approx 10^9$2, and 17–32 pc proper resolution (Ceverino et al., 2018, Ceverino et al., 2021, Nakazato et al., 2024). Snapshot cadence is typically 7–10 Myr in the high-cadence outputs, while the dynamical integration reaches adaptive timesteps down to $M_h \approx 10^9$3 kyr in the sub-$M_h \approx 10^9$4 nebular analyses (Ceverino et al., 2021, Cataldi et al., 6 Oct 2025).

The physics model includes atomic H/He cooling, metal-line cooling, and in several descriptions molecular $M_h \approx 10^9$5 cooling; a cosmological UV background with partial self-shielding; stochastic star formation in gas with $M_h \approx 10^9$6 and $M_h \approx 10^9$7; and multi-channel stellar feedback combining thermal, kinetic, and radiative terms (Ceverino et al., 2017, Ceverino et al., 2018, Cataldi et al., 6 Oct 2025). In the FirstLight I description, radiative feedback is implemented through a non-thermal pressure term in cells containing stars younger than 5 Myr if the gas column density exceeds $M_h \approx 10^9$8, with moderate IR photon trapping above $M_h \approx 10^9$9; momentum injection from supernovae and winds is applied over 40 Myr with IMF-integrated specific momentum $10^{11}\,M_\odot$0 (Ceverino et al., 2017). Metal enrichment from SNe II and Ia is tracked with yields approximating Woosley & Weaver (1995) (Ceverino et al., 2018, Ceverino et al., 2021).

The project’s post-processing stack varies by application. Stellar continua are commonly generated with BPASS; nebular emission is computed either from CLOUDY-based grids or from pixel-by-pixel H II-region models; dust radiative transfer is performed with POLARIS or SKIRT in the dust-focused studies (Ceverino et al., 2018, Ceverino et al., 2021, Mushtaq et al., 2022, Mushtaq et al., 2023, Nakazato et al., 7 Feb 2026). This layered design allows FirstLight to separate hydrodynamical evolution from observable synthesis while retaining the simulated spatial and temporal structure of the ISM.

A recurring limitation is methodological rather than conceptual. H II regions are unresolved in several line-emission papers and are treated through subgrid models with fixed or grid-derived densities; dust studies often adopt a constant dust-to-metal ratio $10^{11}\,M_\odot$1 and fixed MW- or SMC-type grain models; AGN feedback is not included in the standard sub-$10^{11}\,M_\odot$2 analyses and is explicitly stated to be absent in the compactness work (Ceverino et al., 2021, Nakazato et al., 2023, Mushtaq et al., 2023, Cataldi et al., 6 Oct 2025). These choices delimit the regime in which FirstLight predictions should be interpreted.

3. Bursty star formation and redshift-dependent efficiency

A central FirstLight result is that star formation at cosmic dawn is strongly time-variable. In FirstLight II, galaxies spend about 70% of their time in star-formation bursts across $10^{11}\,M_\odot$3 (Ceverino et al., 2018). At $10^{11}\,M_\odot$4, a typical burst has $10^{11}\,M_\odot$5–$10^{11}\,M_\odot$6 and an effective width of about 100 Myr, while roughly one quarter of bursts populate a high-intensity tail with $10^{11}\,M_\odot$7–$10^{11}\,M_\odot$8 and widths of 40–80 Myr (Ceverino et al., 2018). The mean time between consecutive peaks is about $10^{11}\,M_\odot$9 Myr, with only mild mass dependence.

These burst cycles imprint the star-forming main sequence. At $z=5$0, the overall scatter is $z=5$1 dex, rising to $z=5$2 dex at $z=5$3 because of a low-sSFR tail of quiescent dwarfs with $z=5$4 (Ceverino et al., 2018). The mean sSFR grows with redshift approximately as $z=5$5, and the typical burst height remains a factor $z=5$6 above the mean, implying $z=5$7–$z=5$8 at $z=5$9–10, with a tail reaching $z=5$0 (Ceverino et al., 2018).

Later work reframed this burstiness in terms of galaxy formation efficiency. FirstLight defines an instantaneous efficiency

$z=5$1

and an integrated efficiency

$z=5$2

At fixed $z=5$3, the instantaneous efficiency rises from about 10% at $z=5$4 to about 30% at $z=5$5, while the integrated efficiency increases from about 0.1 at $z=5$6 to about 0.2 at $z=5$7 (Ceverino et al., 2024). The same study argues that this redshift dependence is driven by high gas densities in compact early galaxies: efficient systems at $z=5$8 reach $z=5$9, free-fall times $z=5$0 Myr, and $z=5$1 in central regions and off-center clumps (Ceverino et al., 2024).

This combination of high duty cycle, short depletion times, and rising efficiency provides the project’s physical explanation for why galaxies at fixed halo mass are brighter and more rapidly growing at higher redshift. A plausible implication is that FirstLight’s burst statistics are not an incidental by-product of zoom-in numerics, but one of the main mechanisms by which the simulations reproduce the slow evolution of the bright-end UV luminosity function.

4. Spectral energy distributions, metallicity, and nebular diagnostics

FirstLight III established the UV–optical spectral diversity of early galaxies by coupling simulated star-formation histories to BPASS-based stellar populations and CLOUDY nebular grids (Ceverino et al., 2018). At $z=5$2, across $z=5$3, most systems have steep intrinsic UV slopes $z=5$4 to $z=5$5 and high ionizing-photon production efficiencies $z=5$6–25.5, corresponding to $z=5$7–$z=5$8 (Ceverino et al., 2018). Optical lines are correspondingly strong: EW(H$z=5$9+[N II]) commonly lies between 300 and 3000 Å, and broad-band colors in rest-frame $M_{\rm UV} \simeq -12$0, $M_{\rm UV} \simeq -12$1, and $M_{\rm UV} \simeq -12$2 are significantly modified by [O III]+H$M_{\rm UV} \simeq -12$3 and H$M_{\rm UV} \simeq -12$4+[N II] contamination (Ceverino et al., 2018).

The chemical-evolution work predicts that a tight mass–metallicity relation is already present by $M_{\rm UV} \simeq -12$5–8 (Langan et al., 2019). At fixed $M_{\rm UV} \simeq -12$6, the metallicity declines by $M_{\rm UV} \simeq -12$7 dex from $M_{\rm UV} \simeq -12$8 to 5, with a mean decline of about 0.15 dex, similar to the intrinsic scatter (Langan et al., 2019). The scatter correlates with gas fraction: in the $M_{\rm UV} \simeq -12$9 bin, systems with $-22$0 are systematically metal-poor, whereas those with $-22$1 are richer by about 0.2 dex (Langan et al., 2019). The same paper finds that the strong-line index $-22$2 is informative for $-22$3 but saturates at $-22$4 for higher metallicities, and that higher sSFR raises $-22$5 at fixed metallicity through a higher ionization parameter (Langan et al., 2019).

FirstLight IV specialized this picture to sub-$-22$6 galaxies with $-22$7 to $-22$8 at $-22$9–6.5 (Ceverino et al., 2021). In that sample, a factor $z \simeq 5$00 variation in sSFR, from about 1.7 to $z \simeq 5$01, drives a $z \simeq 5$02 dex spread in EW([O III] $z \simeq 5$03): extreme bursts reach EW $z \simeq 5$04 Å, main-sequence systems have EW $z \simeq 5$05 Å, and the most massive, lower-SFR sub-$z \simeq 5$06 galaxies show EW $z \simeq 5$07 Å (Ceverino et al., 2021). Even at fixed sSFR, EW([O III]) retains an intrinsic dispersion $z \simeq 5$08 dex whose residuals correlate with both $z \simeq 5$09 and SFR-weighted $z \simeq 5$10 (Ceverino et al., 2021).

The [O III]-to-H$z \simeq 5$11 ratio is used as a morphological and physical classifier:

$z \simeq 5$12

By convention, [O III]-bright emitters have $z \simeq 5$13, while H$z \simeq 5$14-bright emitters have $z \simeq 5$15 (Ceverino et al., 2021). In the sub-$z \simeq 5$16 sample, [O III]-bright systems tend to have higher ionization parameters and/or higher nebular metallicities, with [O III] up to a factor of $z \simeq 5$17 brighter than H$z \simeq 5$18; H$z \simeq 5$19-bright galaxies dominate the sample and can have H$z \simeq 5$20 up to $z \simeq 5$21 brighter than O III.

The more massive [O III] studies extend the same logic to multi-line diagnostics with ALMA-accessible far-infrared lines (Nakazato et al., 2023). They find that some $z \simeq 5$22 galaxies reach $z \simeq 5$23 already by $z \simeq 5$24, and that high-redshift systems lie above the local $z \simeq 5$25–SFR relation because very young stellar populations boost ionizing-photon production (Nakazato et al., 2023). The line-ratio plane combining [O III] $z \simeq 5$26, [O III] 88 $z \simeq 5$27m, and [O III] 52 $z \simeq 5$28m is then used to infer typical ISM conditions of $z \simeq 5$29–$z \simeq 5$30, $z \simeq 5$31 to $z \simeq 5$32, and $z \simeq 5$33 to $z \simeq 5$34 (Nakazato et al., 2023).

Together, these results make FirstLight a line-physics framework as much as a galaxy-formation suite: stellar age, burst phase, metallicity, and local H II-region conditions are all required to interpret the nebular output of reionization galaxies.

5. Clumps, compactness, and size evolution

The spatial structure of FirstLight galaxies is neither smooth nor static. A merger-driven clump analysis identifies [O III]-bright clumps in 1828 snapshots from 62 galaxies over $z \simeq 5$35–5.5, with clumps found via a threshold $z \simeq 5$36 and minimum radius $z \simeq 5$37 pc (Nakazato et al., 2024). About one tenth of snapshots are clumpy systems with two or more components (Nakazato et al., 2024). Young off-centered clumps have mass-weighted ages of about 35 Myr, radii $z \simeq 5$38 pc, gas fractions $z \simeq 5$39, and sSFR peaks around $z \simeq 5$40, while older clumps in clumpy systems peak around $z \simeq 5$41 and single-component central clumps around $z \simeq 5$42 (Nakazato et al., 2024). The inferred baryonic merger mass ratios are major, $z \simeq 5$43, and the clumpy configurations are short-lived, with off-centered clumps merging into the center on timescales of several tens of Myr, consistent with $z \simeq 5$44 Myr and $z \simeq 5$45 Myr (Nakazato et al., 2024).

A distinct but related line of work studies global compactness. In the compactness analysis, the stellar half-mass radius $z \simeq 5$46 is normalized by the halo radius $z \simeq 5$47 through

$z \simeq 5$48

The distribution of $z \simeq 5$49 is bimodal over $z \simeq 5$50, with median $z \simeq 5$51 and a compact population defined by $z \simeq 5$52 (Cataldi et al., 6 Oct 2025). Galaxies are found to undergo an expansion–compaction–re-expansion cycle. Compaction begins at $z \simeq 5$53 and ends at $z \simeq 5$54, over which the size–mass relation becomes inverted, with a negative slope (Cataldi et al., 6 Oct 2025). The proposed mechanism is a self-reinforced, approximately spherical gas inflow driven by the transition of the central potential from dark-matter dominated to baryon dominated, triggering a strong localized starburst within $z \simeq 5$55 kpc (Cataldi et al., 6 Oct 2025). Unlike low-redshift wet-compaction scenarios, this process ends without central gas depletion or star-formation quenching because the replenishment time becomes comparable to the depletion time (Cataldi et al., 6 Oct 2025).

The later synthetic-observation study extends size evolution to $z \simeq 5$56–5 (Ceverino et al., 5 Mar 2026). In JWST-like rest-frame UV measurements, the size–mass relation is already in place by $z \simeq 5$57 and is fit as

$z \simeq 5$58

For rest-UV sizes, the slope is $z \simeq 5$59 at $z \simeq 5$60–6, $z \simeq 5$61 at $z \simeq 5$62–9, and the normalization evolves rapidly, increasing by about 0.5 dex between $z \simeq 5$63 and $z \simeq 5$64 in 600 Myr (Ceverino et al., 5 Mar 2026). At $z \simeq 5$65, the fitted rest-UV size grows from about 0.17 kpc at $z \simeq 5$66–15 to about 0.60 kpc at $z \simeq 5$67–6 (Ceverino et al., 5 Mar 2026). The same paper finds large diversity at fixed mass: extended galaxies tend to have higher sSFRs and gas fractions and sit above the main sequence, while compact galaxies tend to have lower sSFRs and can occupy “mini-quenching” phases (Ceverino et al., 5 Mar 2026).

The size analysis also shows that dust is dynamically unimportant but observationally decisive. Differential attenuation dims galaxy centers more strongly than outer regions, increasing apparent half-light radii and modifying the size–mass slope even in the rest-frame optical (Ceverino et al., 5 Mar 2026). In the $z \simeq 5$68–6 rest-optical band, the fitted slope changes from $z \simeq 5$69 without dust to $z \simeq 5$70 with dust (Ceverino et al., 5 Mar 2026). This suggests that some of the apparent disagreement between simulations and observed size–mass relations can arise from radiative transfer rather than from dynamics alone.

6. Dust, observability, and the role of FirstLight in JWST-era interpretation

Dust was not included in the intrinsic SED predictions of FirstLight III, but later studies turned it into a major observable channel (Ceverino et al., 2018, Mushtaq et al., 2023). In the POLARIS-based dust-emission study at $z \simeq 5$71 and 8, the CMB sets a floor for diffuse dust temperatures: $z \simeq 5$72 K and $z \simeq 5$73 K (Mushtaq et al., 2022). The resulting equilibrium temperatures exhibit strong spatial gradients, with hot central regions near young stars and outskirts driven toward the CMB floor; including the CMB boosts mid- and far-infrared emission and weakens simple $z \simeq 5$74–$z \simeq 5$75 correlations (Mushtaq et al., 2022).

The attenuation-curve study for $z \simeq 5$76–8 galaxies finds a broad diversity of attenuation laws arising from the clumpy ISM and starbursting geometry (Mushtaq et al., 2023). Low-mass systems with $z \simeq 5$77–$z \simeq 5$78 typically show steep, SMC-like attenuation, whereas more massive galaxies with $z \simeq 5$79 approach Calzetti-like behavior (Mushtaq et al., 2023). At $z \simeq 5$80, the IRX–$z \simeq 5$81 relation follows the Calzetti model with a shift to slightly lower $z \simeq 5$82 at low attenuation due to low metallicity; at $z \simeq 5$83, the locus shifts to higher IRX because of stronger CMB heating (Mushtaq et al., 2023).

A clump-scale extension then decomposes galaxies into clumps, diffuse regions, and system-integrated emission (Nakazato et al., 7 Feb 2026). System-integrated attenuation curves are grayer than Calzetti even for fixed MW- or SMC-type dust, with median slopes $z \simeq 5$84 for MW dust and $z \simeq 5$85 for SMC dust (Nakazato et al., 7 Feb 2026). Individual clumps are even grayer and substantially more obscured: for MW dust, median clump values are $z \simeq 5$86, $z \simeq 5$87 mag, $z \simeq 5$88 mag, and $z \simeq 5$89 (Nakazato et al., 7 Feb 2026). Diffuse regions are steeper and optically thinner, with $z \simeq 5$90, $z \simeq 5$91 mag, and $z \simeq 5$92 (Nakazato et al., 7 Feb 2026). The paper interprets this through a toy model on the IRX–$z \simeq 5$93 plane, parameterized by the dust-to-star scale-height ratio $z \simeq 5$94 and a fiducial optical depth $z \simeq 5$95 (Nakazato et al., 7 Feb 2026).

These radiative-transfer outputs feed directly into instrumental predictions. For the sub-$z \simeq 5$96 sample at $z \simeq 5$97, galaxies with $z \simeq 5$98 have observed $z \simeq 5$99-band magnitudes 27–28 and typical optical-line fluxes of a few $h^{-1}\,\mathrm{Mpc}$00, making them accessible to JWST/NIRSpec in 0.5–2 hour exposures (Ceverino et al., 2021). FirstLight IV predicts that [O III]-bright emitters will often be barely resolved at JWST’s nominal FWHM $h^{-1}\,\mathrm{Mpc}$01 arcsec, whereas H$h^{-1}\,\mathrm{Mpc}$02 emission should be better resolved and more extended where it is observable (Ceverino et al., 2021).

The MIRI-oriented spectroscopy study uses the FirstLight database to forecast detectability at $h^{-1}\,\mathrm{Mpc}$03 (Álvarez-Márquez et al., 2019). In 40 ks with MIRI/MRS, only the luminous tail is detectable at $h^{-1}\,\mathrm{Mpc}$04: 6.2%, 1.1%, and 0.4% of the simulated population in the redshift bins 6.5–7.5, 7.5–8.5, and $h^{-1}\,\mathrm{Mpc}$05, corresponding to intrinsic limits of $h^{-1}\,\mathrm{Mpc}$06, 1.9, and $h^{-1}\,\mathrm{Mpc}$07 and stellar masses $h^{-1}\,\mathrm{Mpc}$08, 9, and $h^{-1}\,\mathrm{Mpc}$09, respectively (Álvarez-Márquez et al., 2019). For bright or lensed sources, the same study finds that 10–40 ks spectra can recover H$h^{-1}\,\mathrm{Mpc}$10, [O III] $h^{-1}\,\mathrm{Mpc}$11, H$h^{-1}\,\mathrm{Mpc}$12, He I 1.083 $h^{-1}\,\mathrm{Mpc}$13m, and, in deeper integrations, [N II], [S II], and [S III], enabling extinction, excitation, hardness, and metallicity diagnostics (Álvarez-Márquez et al., 2019).

The broader significance of FirstLight lies in this observational closure. The suite was built to reproduce the UV luminosity function and scaling relations of early galaxies, but it now functions as a multi-layer inference engine linking bursty accretion, feedback-regulated star formation, chemical enrichment, nebular line physics, dust geometry, compactness, and synthetic observables. This suggests that its enduring value is not any single prediction, but the internally connected description it provides of how reionization-era galaxies grow, radiate, and appear in data.