1. Introduction
Recent estimates of the star formation histories of high-redshift galaxies using the James Webb Space Telescope (JWST) indicate that star formation proceeds in bursts (i.e. is rapidly varying) in the early Universe (e.g. Faisst et al. Reference Faisst, Capak, Emami, Tacchella and Larson2019; Looser et al. Reference Looser2023; Strait et al. Reference Strait2023). Bursty star formation occurs when the effective equilibrium between the self-gravity of gas in a galaxy and stellar feedback following star formation cannot be maintained (cf. Faucher-Giguère Reference Faucher-Giguère2018). It is expected to be commonplace in lower-mass galaxies (e.g. Oñorbe et al. Reference Oñorbe2015; Muratov et al. Reference Muratov2015; Sparre et al. Reference Sparre2017) and, importantly, in galaxies at high redshifts (e.g. Pallottini & Ferrara Reference Pallottini and Ferrara2023; Shen et al. Reference Shen2023; Sun et al. Reference Sun, Faucher-Giguère, Hayward and Shen2023).
Previous work has shown that bursty star formation will produce naturally strong stellar feedback-driven outflows in high-z galaxies, especially at lower masses (e.g. Furlanetto & Mirocha Reference Furlanetto and Mirocha2022; Menon & Power Reference Menon and Power2024). This results in periods of oscillating gas mass in the galaxy as gas is expelled from the interstellar medium (ISM) and cosmological gas accretion from the intergalactic medium (IGM) is suppressed. The massive stars that are the sources of feedback, via e.g. stellar winds and supernovae (e.g. Lamers & Cassinelli 1999), are also the sources of metals that enrich the ISM (e.g. Mannucci et al. Reference Mannucci, Cresci, Maiolino, Marconi and Gnerucci2010; Maiolino & Mannucci Reference Maiolino and Mannucci2019) as well as the primary dust production mechanism in the early Universe (e.g. Todini & Ferrara Reference Todini and Ferrara2001; Gall & Hjorth Reference Gall and Hjorth2018).
This prompts the important question of how bursty star formation influences the metal and dust content of high-redshift galaxies. From a theoretical perspective, we expect that metals and dust will play a crucial role in galaxy formation as coolants (e.g. Cox & Tucker Reference Cox and Tucker1969; Sutherland & Dopita Reference Sutherland and Dopita1993; Ploeckinger & Schaye Reference Ploeckinger and Schaye2020), which plays a critical role in regulating a galaxy’s star formation and associated feedback efficiencies. From an observational perspective, a galaxy’s metallicity will influence what we can infer from its spectral energy distribution (e.g. Robotham et al. Reference Robotham2020; Vijayan et al. Reference Vijayan2025), including diagnostics such as star formation rate. The presence of dust will obscure star forming regions in individual galaxies and star forming galaxies (e.g. Adelberger & Steidel Reference Adelberger and Steidel2000); there are indications from ALMA (Atacama Large Millimeter/submillimeter Array) and JWST data that galaxies at
$z\geq\,6$
harbour a diverse range of dust content (e.g. Matthee et al. Reference Matthee2019; Rodighiero et al. Reference Rodighiero2023; Barrufet et al. Reference Barrufet2023).
Astrophysically, the heavy elements that drive the growth of metallicity are the product of nucleosynthesis in the cores of stars, supernovae, and stellar mergers (cf. Nomoto, Kobayashi, & Tominaga Reference Nomoto, Kobayashi and Tominaga2013; Kobayashi, Karakas, & Lugaro Reference Kobayashi, Karakas and Lugaro2020; Arcones & Thielemann Reference Arcones and Thielemann2022), which subsequently enrich the ISM, circumgalactic medium (CGM), and IGM via stellar-driven winds and supernovae. The physics of dust formation and evolution is complex (e.g. McKee Reference McKee1989; Dwek Reference Dwek1998; Calura Reference Calura2025), but its production is believed to be driven primarily by supernovae at redshifts
$z\geq5$
(e.g. Todini & Ferrara Reference Todini and Ferrara2001; Gall, Hjorth, & Andersen Reference Gall, Hjorth and Andersen2011; Gall & Hjorth Reference Gall and Hjorth2018; Leśniewska & Michał owski Reference Leśniewska and Michałowski2019; Ferrara et al. Reference Ferrara2022; Langeroodi et al. Reference Langeroodi, Hjorth, Ferrara and Gall2024). Feedback that drives gas out of a galaxy will lower the gas-phase metallicity of the ISM, whereas fresh accretion of nearly pristine gas from the IGM will lower the ISM gas phase metallicity. Similarly, supernovae-driven shocks in the ISM destroy dust, while they will also be entrained in the metal-enriched gas that is expelled from galaxies via winds. As a result, we might expect to see differences with respect to star formation modes in which burstiness is suppressed. This is because galaxies can retain more gas and more metal-enriched gas and dust in their ISM, which should allow for the gas-metallicity and dust content to increase steadily with time. Recent work by Marszewski et al. (Reference Marszewski, Faucher-Giguère, Feldmann and Sun2025) using the FIRE-2 simulations suggests that bursty star formation is necessary to reproduce the nearly constant mass-metallicity relation at
$z\geq5$
– strong outflows induced by bursty star formation act to reset the ISM in a galaxy, and so metal-enriched accretion onto the galaxy is offset by reduced metal-production efficiency in the ISM. As Liu et al. (Reference Liu, Sun, Chang, Furlanetto and Bradford2024) show, bursty star formation’s impact on metallicity, predominantly in lower-mass galaxies, introduces a strongly mass-dependent scatter that helps to explain ALMA [CII] metal line intensity mapping measurements during the Epoch of Reionisation.
Interestingly, we might expect the efficiency of feedback to be influenced by metallicity at early epochs. The first generation of stars, the zero-metallicity Population III, enriched the ISM of their host galaxies (Sluder et al. Reference Sluder, Ritter, Safranek-Shrader, Milosavljević and Bromm2015; Chen et al. Reference Chen2024) as supernovae (Heger & Woosley Reference Heger and Woosley2002, Reference Heger and Woosley2010), which led to subsequent generations of Population II stars being metal enriched. We expect increased stellar metallicity to influence stellar mass loss rates via winds, and consequently the efficiency of stellar feedback (e.g. Krtička & Kubát 2014; Gormaz-Matamala et al. Reference Gormaz-Matamala2022; Rickard et al. Reference Rickard2022; Dekel et al. Reference Dekel, Sarkar, Birnboim, Mandelker and Li2023). It also regulates the maximum progenitor mass above which stars collapse directly into black holes without exploding as supernova (e.g. Zhang, Woosley, & Heger Reference Zhang, Woosley and Heger2008; O’Connor & Ott Reference O’Connor and Ott2011; Jecmen & Oey Reference Jecmen and Oey2023), which reduces the stellar feedback efficiency for a given stellar population.
In this paper, we study how bursty star formation influences the metal and dust content of high-z galaxies - at
$z\geq 5$
, corresponding to the Cosmic Dawn and the Epoch of Reionisation – using our lightweightFootnote
a
semi-analytical model Ashvini
Footnote
b
(see Menon & Power Reference Menon and Power2024, hereafter Reference Menon and PowerMP24). We extend Ashvini to track the growth of gas-phase and stellar metallicities, and dust, and we also include a model to account for the metallicity-dependence of the efficiency of feedback. These additions allow us to explore how bursty star formation during reionisation imprints scatter and asymmetry in mass–metallicity relations, and how feedback cycles delay chemical and dust enrichment in halos of different mass ranges. We do this by investigating how our results are sensitive to key model parameters, and comparing to a model in which feedback is instantaneous and star formation proceeds in a smooth (i.e. non-bursty) manner.
The paper is structured as follows. In Section 2, we briefly review Ashvini as implemented in Reference Menon and PowerMP24 before providing a detailed description on our updates to metal (Section 2.2.2) and dust (Section 2.2.4) evolution, and metallicity-dependent feedback (Section 2.2.3). In Section 3, we present our main results, reviewing our predictions for gas-phase and stellar metallicity growth over cosmic time; the relationship between metallicity and stellar mass; and predictions for dust growth over time. We review these results in the context of previous observational and theoretical work in Section 4 and summarise our conclusions in Section 5.
Note that we use the following values for the cosmological parameters:
$\Omega_{\textrm{b}}=0.0484$
,
$\Omega_{\textrm{M}}=0.308$
,
$\Omega_{\Lambda}=0.692$
,
$h=0.678$
,
$\sigma_8=0.815$
, and
$n_{\textrm{s}}=0.968$
, which are consistent with the results obtained by Planck Collaboration et al. (2020). We assume a metal mass fraction in the solar neighbourhood of
${Z}_{\odot}=0.015$
(Asplund et al. Reference Asplund, Grevesse, Sauval and Scott2009; Lodders Reference Lodders2019).
2. Theoretical model
We use our lightweight semi-analytical model Ashvini, Footnote c introduced in Reference Menon and PowerMP24, to study the evolution of gas and stellar components in dark matter halos, which is motivated by the equilibrium model approach of Davé, Finlator, & Oppenheimer (Reference Davé, Finlator and Oppenheimer2012). The model allows for delays between when stars form and when they produce feedback (e.g. Furlanetto & Mirocha Reference Furlanetto and Mirocha2022) and time-dependent suppression of cosmological gas accretion, which reflects the growth of an ionising UV background (UVB) in the early Universe (e.g. Kravtsov & Manwadkar Reference Kravtsov and Manwadkar2022). Note that we do not account for the effects of the high-z progenitors of super-massive black holes and the feedback they produce; this will be incorporated in a forthcoming update to Ashvini.
2.1. Input halo merger trees
We use the same merger trees as in Reference Menon and PowerMP24, which were generated with the Monte Carlo algorithm of Parkinson, Cole, & Helly (Reference Parkinson, Cole and Helly2008). Each merger tree consists of a halo mass assembly history in the redshift range
$5\leq z\leq 25$
, equally spaced in the logarithm of the expansion factor,
$a = 1/(1 + z)$
. We consider 11 mass bins at
$z=5$
that are spaced in a quasi-logarithmic progression, alternating between factors of 10 and 5, within the mass range
$10^6\leq M_{\textrm{h}}/\text{M}_{\odot} \leq 10^{11}$
. We generate merger trees for 100 halos in each mass bin, such that we have 100 halos with masses of e.g.
$M_{\textrm{h}}=10^6\,\text{M}_{\odot}$
at
$z=5$
. Each merger tree provides us with
$M_{\textrm{h}}$
as a function of z, and so we can calculate
$\dot{M}_{\textrm{h}}$
as a function of redshift (z) or time (
$t\equiv t(z)$
) by constructing a smooth spline interpolant, which we differentiate numerically. In Reference Menon and PowerMP24, we had shown that the mass bin
$10^7\,\text{M}_{\odot}$
marks the transition from ‘low’ to ‘high’ mass halos. Low/high-mass halos are particularly susceptible/robust to feedback processes due to their shallow/deep potential wells.
2.2. The Ashvini model
In this subsection, we expand on the new updates to the Ashvini model. First, we give a brief summary of the model; we refer the reader to Reference Menon and PowerMP24 for a more detailed discussion.
2.2.1. Core model
Ashvini employs a set of coupled differential equations to describe the evolution of gas and stellar mass within a dark matter halo across cosmic time, and enforces mass conservation by balancing inflows, outflows, and internal processes. At a given instant in time, t, we compute the growth rate of the gas mass
$\dot{M}_{\textrm{g}}$
, the star formation rate
$\dot{M}_{\star}$
, and the wind-driven mass loss rate
$\dot{M}_{\textrm{w}}$
from the following:
\begin{align}\dot{M}_{\textrm{g}}&=\dot{M}_{\textrm{c},\textrm{g}}-\dot{M}_{\star} - \dot{M}_{\textrm{w}}(t^\prime),\\[-10pt] \nonumber \end{align}
\begin{align}\dot{M}_{\star}&=\epsilon_{\text{sf}}\frac{M_{\textrm{g}}}{\tau_{\text{sf}}}, \\[-10pt] \nonumber \end{align}
\begin{align}\dot{M}_{\textrm{w}}&=\eta_{\text{fb}}\dot{M}_{\star}(t^\prime), \\[12pt] \nonumber \end{align}
where we assume the
$\dot{M}$
values correspond to the instantaneous rates at t.
$t^\prime = t - t^{\textrm{d}}$
allows for a time delay (
$t^{\textrm{d}}\gt0$
) between when stars form and when they produce feedback;
$t^{\textrm{d}}=0$
corresponds to the case of instantaneous feedback. In Equation (1),
$\dot{M}_{\textrm{c},\textrm{g}}$
is the cosmological gas accretion from the IGM (see below); in Equation (2),
$M_{\textrm{g}}$
is the instantaneous mass of the gas reservoir, and
$\epsilon_{\text{sf}}$
and
$\tau_{\text{sf}}$
are the star formation efficiency and the star formation timescale (see below) respectively; and in Equation (3),
$\eta_{\text{fb}}$
is the feedback efficiency.
As in Reference Menon and PowerMP24, we estimate the cosmological gas accretion rate using,
\begin{equation*}\dot{M}_{\textrm{c},\textrm{g}}=\varepsilon_{\mathrm{UV}}\bigg(\frac{\Omega_{\textrm{b}}}{\Omega_{\textrm{M}}}\bigg)\dot{M}_{\textrm{h}},\end{equation*}
where
$\dot{M}_{\textrm{h}}$
is the halo mass growth rate;
$(\Omega_{\textrm{b}}/\Omega_{\textrm{M}})$
is the cosmic baryon fraction; and
$0\leq\varepsilon_{\text{UV}}\leq1$
quantifies the degree to which the accretion rate is suppressed by the UVB, which we assume is effective from a reionisation redshift,
$z_{\text{rei}}{({=}7)}$
. The specific value for
$\varepsilon_{\text{UV}}$
at a given z is fixed by our adopted model for UVB suppression from Reference Menon and PowerMP24, as summarised in Appendix 1. We also assume that
$\tau_{\text{sf}}$
is a fraction of the Hubble time,
\begin{equation*} \tau_{\text{sf}} = 0.15 \frac{f_{\text{sf}}}{H(z)}.\end{equation*}
We set
$f_{\text{sf}}$
= 1 for simplicity. In practice, we expect
$f_{\text{sf}}\lt1$
in local star forming regions where dynamical times will be much shorter than the halo-averaged value, but
$f_{\text{sf}}$
provides a reasonable upper limit given the uncertainties.
2.2.2. Modelling gas and star metallicity
In this work, we update the Ashvini model to track the growth of the mass of metals in gas and stars,
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
. We extend Equations (1)–(3) as:
\begin{align}\dot{M}_{\textrm{Z},\textrm{g}}&=Z_{\text{IGM}}\dot{M}_{\textrm{c},\textrm{g}} + Y_{\textrm{Z}}\dot{M}_{\star}(t^\prime)-\dot{M}_{\textrm{Z},\star}-\dot{M}_{\textrm{Z},\textrm{w}}(t^\prime), \\[-10pt] \nonumber \end{align}
\begin{align}\dot{M}_{\textrm{Z},\star}&=\epsilon_{\text{sf}}\dfrac{M_{\textrm{Z},\textrm{g}}}{\tau_{\text{sf}}}, \\[-10pt] \nonumber\end{align}
\begin{align}\dot{M}_{\textrm{Z},\textrm{w}}&=\frac{M_{\textrm{Z},\textrm{g}}}{M_{\textrm{g}}}\dot{M}_{\textrm{w}}, \\[12pt] \nonumber \end{align}
where
$Z_{\text{IGM}}$
is the metallicity of the freshly accreted IGM gas and
$Y_{\textrm{Z}}$
is the heavy element yield per unit star formation rate (SFR).
-
• Equation (4) tracks the metal content of the accreting gas, which we assume to be at a fixed value
$Z_{\text{IGM}}$
times the cosmological mass accretion rate
$\dot{M}_{\textrm{c},\textrm{g}}$
; the metals formed in the cores of stars that would eventually enrich the ISM once the star dies,
$Y_{\textrm{Z}}\dot{M}_{\star}(t-t_{\textrm{d}})$
; the metals in the ISM that are locked in newly formed stars,
$\dot{M}_{\textrm{Z},\star}$
; and the metals that are lost via winds
$\dot{M}_{\textrm{Z},\textrm{w}}$
. -
• Equation (5) assumes that the mass of metals in stars tracks the star formation rate
$\dot{M}_{\ast}$
as given by Equation (2), times
$(M_{\textrm{Z},\textrm{g}}/M_{\textrm{g}})$
, at that instant of cosmic time. -
• Equation (6) assumes that the gas mass of metals lost via winds is simply
$(M_{\textrm{Z},\textrm{g}}/M_{\textrm{g}})$
times
$\dot{M}_{\textrm{w}}$
at that instant, which is driven by the star formation rate
$\dot{M}_{\ast}$
at the earlier time,
$t^\prime$
.
These equations assume an instantaneous perfect mixing of metals with the ambient gas reservoir. This is a simplifying assumption, but it is consistent with the rapid mixing timescales expected in a turbulent, feedback-driven ISM (Pan, Scannapieco, & Scalo Reference Pan, Scannapieco and Scalo2013; Hirai & Saitoh Reference Hirai and Saitoh2017) that we would expect in high redshift galaxies.
We use as our fiducial values
$Z_{\text{IGM}}=10^{-3}\,\text{Z}_{\odot}$
, following Kravtsov & Manwadkar (Reference Kravtsov and Manwadkar2022), and
$Y_{\textrm{Z}}=0.06$
, following Vincenzo et al. (Reference Vincenzo, Matteucci, Belfiore and Maiolino2016) assuming a Chabrier (Reference Chabrier2003) initial mass function (IMF).
2.2.3. Modelling metal-dependent feedback
We employ a momentum-regulated feedback recipe that balances the momentum injected into the ISM gas, by supernovae as well as other feedback processes, to that which is required to eject the gas at the halo’s escape velocity (also see Furlanetto et al. Reference Furlanetto, Mirocha, Mebane and Sun2017; Furlanetto & Mirocha Reference Furlanetto and Mirocha2022). This yields a prescription that depends on the redshift, z, and the halo mass,
$M_{\textrm{h}}$
as,
\begin{equation}\eta_{\text{fb}}=\epsilon_{\text{fb}}\pi_{\textrm{p}}f(Z_{\star})\bigg(\frac{10^{11.5}\,\text{M}_{\odot}}{M_{\textrm{h}}}\bigg)^{1/3}\bigg(\frac{9}{1+z}\bigg)^{1/2}.\end{equation}
Here
$\pi_{\textrm{p}}$
(which we set to unity) is the total momentum per supernova (in units of
$2\times10^{33}$
g cm s
$^{-2}$
, see also Furlanetto et al. Reference Furlanetto, Mirocha, Mebane and Sun2017; Furlanetto & Mirocha Reference Furlanetto and Mirocha2022; Menon & Power Reference Menon and Power2024) and
$\epsilon_{\mathrm{fb}}$
is the fraction of this momentum that is coupled to the galaxy’s ISM as a wind.
The factor
$f(Z_{\ast})\leq1$
governs how the feedback efficiency depends on the stellar metallicity,
$Z_{\ast}$
. If we assume that there is no metal dependence, then
$f(Z_{\ast})=1$
. However, we expect that the feedback should be metal-dependent at low stellar metallicities (eg. Zhang et al. Reference Zhang, Woosley and Heger2008; O’Connor & Ott Reference O’Connor and Ott2011; Sukhbold et al. Reference Sukhbold, Ertl, Woosley, Brown and Janka2016; Jecmen & Oey Reference Jecmen and Oey2023). In our model, we adopt an empirical parameterisation for
$f(Z_{\star})$
motivated by the work of Jecmen & Oey (Reference Jecmen and Oey2023). They estimated reductions in the total integrated momentum and mechanical energy of 75% and 40% respectively for stellar metallicities
$Z_{\star}\lesssim 0.4\,\text{Z}_{\odot}$
compared to values expected at solar metallicity. We model this effect by means of a sigmoid function, which produces step-like behaviour with a smooth transition, given by,
\begin{equation} f(Z_{\star})=\frac{\left(b-a\right)}{1+\exp-((Z_{\star}-m)/s)}+a. \end{equation}
Here,
$m=0.4\,\text{Z}_{\odot},s=0.0067\,\text{Z}_{\odot},a=0.25$
and
$b=1$
are parameters that control the amplitude and sharpness of the transition of the sigmoid; this gives
$f(Z_{\star})=0.25$
at zero-metallicity.
2.2.4. Modelling dust
Although the physics of dust formation and evolution is complex (e.g. McKee Reference McKee1989; Dwek Reference Dwek1998; Calura Reference Calura2025), we incorporate an approximate model to obtain estimates of the dust mass associated with the gas phase,
$M_{\textrm{d}}$
. At redshifts
$z\geq5$
, dust production is believed to be driven primarily by supernovae (e.g. Todini & Ferrara Reference Todini and Ferrara2001; Gall et al. Reference Gall, Hjorth and Andersen2011; Gall & Hjorth Reference Gall and Hjorth2018; Leśniewska & Michał owski 2019; Ferrara et al. Reference Ferrara2022; Langeroodi et al. Reference Langeroodi, Hjorth, Ferrara and Gall2024) because channels linked to evolved lower mass stars (e.g. the asymptotic giant branch, AGB) are unimportant until later cosmic times (Marassi et al. Reference Marassi2019; Triani et al. Reference Triani, Sinha, Croton, Pacifici and Dwek2020). The standard assumption is that dust is produced by core-collapse supernovae of stars in the mass range 8–40 M
$_{\odot}$
(cf. Heger et al. Reference Heger, Fryer, Woosley, Langer and Hartmann2003; Gall & Hjorth Reference Gall and Hjorth2018) and the dust yield per supernova is a function of the metallicity and mass of the progenitor mass. Gall & Hjorth (Reference Gall and Hjorth2018) estimate this dust yield per supernova as
$0.31\,\text{M}_{\odot}$
with an uncertainty of 0.15 dex. Higher redshift galaxies experience more efficient dust removal than their lower redshift counterparts, driven by active galactic nucleus (AGN) activity, supernova shocks, and astration, with recent estimates of dust removal rates showing a decline from 1.8 Gyrs at
$z\simeq0.05$
to 450 Myrs at
$z\geq3$
(cf. Leśniewska, Hjorth, & Gall Reference Leśniewska, Hjorth and Gall2025).
In Ashvini, we assume that the rate of dust mass growth tracks the rate of core-collapse supernovae, which is proportional to
$\dot{M}_{\star}$
, and is lost through destruction in supernova shocks
$(\dot{M}_{\textrm{d},\mathrm{dest}})$
and expulsion via winds
$(\dot{M}_{\textrm{w}})$
. We describe this using the equations,
\begin{equation}\dot{M}_{\textrm{d}}=Y_{\textrm{d}}\dot{M}_{\star}(t^\prime)- \dot{M}_{\textrm{d}, \mathrm{ dest}}(t^\prime)-\dot{M}_{\textrm{d},\textrm{w}}(t^\prime),\end{equation}
with
\begin{equation}\dot{M}_{\textrm{d},\textrm{w}}=\left(\frac{M_{\textrm{d}}}{M_{\textrm{g}}}\right)\dot{M}_{\textrm{w}}.\end{equation}
Here, the total dust yield per unit star formation rate is
$Y_{\textrm{d}}=0.004$
for our assumed Chabrier (Reference Chabrier2003) IMF, following Gall & Hjorth (Reference Gall and Hjorth2018) and Langeroodi et al. (Reference Langeroodi, Hjorth, Ferrara and Gall2024). The prefactor
$(M_{\textrm{d}}/M_{\textrm{g}})$
in Equation (10) assumes that dust mass loss via winds is proportional to the gas mass loss via winds – we assume an instantaneous perfect mixing of dust with the ambient gas reservoir – and
$\dot{M}_{\textrm{d}}, \mathrm{dest}$
determines the rate at which dust is destroyed by supernovae shocks.
The most straightforward parameterisation for
$\dot{M}_{\textrm{d}, \mathrm{dest}}$
is,
\begin{equation}\dot{M}_{\textrm{d}, \mathrm{dest}}=\frac{M_{\textrm{d}}}{\tau_{\text{dest}}},\end{equation}
where
$\tau_{\text{dest}}$
is the dust destruction timescale (e.g. Triani et al. Reference Triani, Sinha, Croton, Pacifici and Dwek2020). While Equation (11) is the conventional form, we instead adopt a parameterisation that explicitly links to the physical quantities important for dust destruction in our model. Because we assume that supernovae shocks drive dust destruction, we make explicit the dependence on the supernovae rate, which will govern the dust destruction rate, and choose,
\begin{equation} \dot{M}_{\textrm{d}, \mathrm{dest}}= \epsilon_{\textrm{d}}\,M_{\text{SN,sw}}R_{\text{SN}}\,f(M_{\textrm{g}}/M_{\text{crit}}). \end{equation}
Here
$R_{\text{SN}}$
is the rate of supernovae,
$M_{\text{SN,sw}}$
is the mass of the ISM swept up by supernovae ejecta, and
$0\leq\epsilon_{\textrm{d}}\leq1$
governs the mass of dust in the swept-up ISM that is destroyed. In practice, we fix the product
$\epsilon_{\textrm{d}} M_{\text{SN,sw}}$
, as we show below. The function
$f(M_{\textrm{g}}/M_{\text{crit}})$
allows for supernovae to either sweep out dust from the potential in lower-mass galaxies or destroy dust in the ISM in higher-mass galaxies;
$M_{\text{crit}}$
controls the transition between lower- and higher-mass galaxies.
Estimating
$R_{\textrm{SN}}$
: We evaluate
$R_{\text{SN}}$
in our model at time t using the star formation rate (
$\dot{M}_{\star}$
) and the number of core collapse supernovae per unit stellar mass formed,
$\gamma$
.
$\dot{M}_{\star}$
is estimated at
$t^\prime$
to account for the delay between when stars form and when they produce feedback.
$\gamma$
can be calculated directly from the chosen IMF (
$\xi$
) as,
\begin{equation} \gamma = \frac{\int_8^{40}\xi(m)dm}{\int_{0.1}^{100}m\xi(m)dm}.\end{equation}
Taken together, this gives,
\begin{equation}R_{\text{SN}}=\gamma{\dot{M}_{\star}(t^\prime)},\end{equation}
where
$\gamma$
varies between
$1.014\times10^{-2}\,\text{M}_{\odot}^{-1}$
for a Kroupa (Reference Kroupa2001) IMF, to
$1.077\times10^{-2}\,\text{M}_{\odot}^{-1}$
for a Chabrier (Reference Chabrier2003) IMF, to
$1.958\times10^{-2}\,\text{M}_{\odot}^{-1}$
for a Salpeter (Reference Salpeter1955) IMF. We assume the Chabrier (Reference Chabrier2003) as our default value in this paper.
Estimating
${\boldsymbol{\epsilon}}_{{\textbf{d}}}{\textbf{{M}}}_{\textbf{{SN,sw}}}$
: Previous work has estimated that
$M_{\text{SN,sw}}$
is of order the amount of swept-up mass in the warm neutral medium (WNM) during the initial phase of a supernova, when the ejecta velocity is
$\sim$
$10^4$
km s
$^{-1}$
, until it enters the Sedov phase. The WNM has temperatures and number densities
$T\simeq10^4\,\text{K}$
and
$n\sim0.1\,\text{cm}^{-3}$
, and dust destruction occurs as a result of collisions with high-velocity ions (
$\geq$
100 km s
$^{-1}$
). For a multiphase ISM,
$M_{\text{SN,sw}}\lesssim10^4\,\text{M}_{\odot}$
while
$\epsilon_{\textrm{d}}=0.2$
(e.g. McKee Reference McKee1989; Calura Reference Calura2025). For simplicity, we assume that
$\epsilon_{\textrm{d}}\,M_{\text{SN,sw}}=10^3\,\text{M}_{\odot}$
(e.g. Gall & Hjorth Reference Gall and Hjorth2018; Calura Reference Calura2025); the precise value affects the normalisation but not the qualitative behaviour of dust evolution.
Estimating
$\textbf{{f}}({\textbf{{M}}}_{\textbf{{g}}}/{\textbf{{M}}}_{\textbf{{crit}}})$
: Dust destruction depends on whether the galaxy retains sufficient gas mass for supernovae to efficiently couple their energy to the ISM. We use the total gas mass
$M_{\textrm{g}}$
as a threshold: in low-mass systems (
$M_{\textrm{g}} \lt M_{\text{crit}}\simeq10^5\,\text{M}_{\odot}$
), supernova remnants break out of the galaxy before entering the Sedov-Taylor phase, ejecting dust via winds rather than destroying it in situ. We set
$f(M_{\textrm{g}}/M_{\text{crit}})\rightarrow0$
in this regime. In high-mass systems (
$M_{\textrm{g}}\geq M_{\text{crit}}$
), remnants sweep up significant mass, efficiently destroying dust via shocks, and
$f(M_{\textrm{g}}/M_{\text{crit}})\rightarrow1$
. To capture the smooth transition between these regimes, we adopt:
\begin{equation}f(M_{\textrm{g}}/M_{\text{crit}}) = 1 - e^{-\,\left(M_{\textrm{g}}/{M_{\text{crit}}}\right)^\alpha}, \end{equation}
with
$\alpha=8$
producing a sharp transition over
$\sim$
1 dex in mass. This allows for a range over which supernovae effects on dust shift from ejective to destructive. Dust mass in our model grows in proportion to the star formation rate (via core collapse supernovae), is destroyed by supernovae shocks in sufficiently massive galaxies, and is lost via winds in low-mass systems.
3. Results
In this section, we present results for the trends in (gas and stellar) metallicities and dust mass with the Ashvini model. We also explore the impact of different free parameters in our models on these trends. Table 1 lists our model parameters along with their fiducial values.
Parameter values for our fiducial model.
3.1. Evolution of
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
In this subsection, we focus on model predictions for the mass bin
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
at
$z=5$
. We do this by tracking the behaviour of the median and
$10{\text{th}}$
-to-
$90{\text{th}}$
percentile variation for gas and stellar mass,
$M_{\textrm{g}}$
and
$M_{\star}$
, and mass of metals in gas and stars,
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
, as a function of cosmic time, using the Monte Carlo trees described in Section 2.1. We assume the fiducial model parameters as given in Table 1. In those cases in which we show model predictions based on parameter variations (e.g. IGM metallicity,
$Z_{\text{IGM}}$
), we show the predictions of the fiducial model as light-grey curves. As argued in Reference Menon and PowerMP24,
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
is an interesting mass bin to examine because it marks the transition between low-mass halos that are quenched by delayed feedback alone and high-mass halos that continue to accrete gas and form stars to later times.
Fiducial Model Predictions: We begin by assessing how the growth of
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
is influenced by delayed feedback and the resulting bursty star formation. We compare our model predictions for the cases of instantaneous and delayed feedback, with the UV suppression of gas accretion in both cases.
Figure 1 shows our model predictions assuming fiducial parameters for the instantaneous (delayed) feedback scenario (left and right panels) for halo mass bins of
$M_{\textrm{h}} = 10^7\,\text{M}_{\odot}$
and
$10^{10}\,\text{M}_{\odot}$
(upper and lower panels). The solid and dot-dashed curves correspond to the median values of
$M_{\textrm{g}}$
and
$M_{\star}$
; dashed and dotted curves correspond to the median values of
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
; and the shaded regions indicate the
$10{\text{th}}$
-to-
$90{\text{th}}$
percentile variations in each quantity.
Impact of instantaneous versus delayed feedback (fiducial model parameters): The evolution of gas and stellar mass (
$M_{\textrm{g}}$
and
$M_{\star}$
; solid and dotted-dashed curves), and the mass of metals in the gas phase and stars (
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
; dashed and dotted curves) as a function of cosmic time (in Gyrs, lower horizontal axis) and redshift (upper horizontal axis). These predictions are based on the assembly histories of 100 halos with
$M_{\textrm{h}} = 10^7\,\text{M}_{\odot}$
(
$M_{\textrm{h}} = 10^{10}\,\text{M}_{\odot}$
) mass bin at
$z=5$
in the upper (lower) panel, for instantaneous (delayed) feedback in the left (right) panel. Each curve represents the evolution of the median value of a given quantity, while bands indicate the range of the 10th and 90th percentiles.
We first focus on the behaviour in the
$M_{\textrm{h}} = 10^7\,\text{M}_{\odot}$
mass bin (upper panels). The case of instantaneous feedback shows that the evolution of
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
closely follows that of
$M_{\textrm{g}}$
and
$M_{\star}$
, respectively, with similar shapes but offset in normalisation. Both
$M_{\textrm{g}}$
and
$M_{\textrm{Z},\textrm{g}}$
show a smooth monotonic increase to their maximum value at
$z\simeq7$
, which corresponds to the redshift of reionisation,
$z_{\text{rei}}$
, before a subsequent decline of
$\sim$
3 dex by
$z\simeq 5$
.
$M_{\star}$
and
$M_{\textrm{Z},{\star}}$
show a steady increase with cosmic time. The ratio of the two sets of curves remains consistent at
$M_{\textrm{Z},\textrm{g}}/M_{\textrm{g}}=M_{\textrm{Z},{\star}}/{M}_{\star}\simeq10^{-4}$
.
The case of delayed feedback shows some interesting differences with respect to the trends in the instantaneous feedback case. Both
$M_{\textrm{g}}$
and
$M_{\textrm{Z},\textrm{g}}$
show oscillatory behaviour at early times that closely track each other, although the oscillations are more pronounced in
$M_{\textrm{g}}$
. The median
$M_{\textrm{g}}$
and
$M_{\textrm{Z},\textrm{g}}$
fall sharply to zero at
$z\simeq 6.5$
, although the
$10{\text{th}}$
-to-
$90{\text{th}}$
percentile variations reveal that the metal-enriched gas persists in a subset of the halo sample to
$z\simeq5.5$
. There are imprints of the oscillations in
$M_{\textrm{g}}$
evident in
$M_{\star}$
and
$M_{\textrm{Z},{\star}}$
at early times, but otherwise they show a steady increase in cosmic time. Despite these differences, we see that the ratio of the two sets of curves mirrors that in the instantaneous case and remains consistent at
$M_{\textrm{Z},\textrm{g}}/M_{\textrm{g}}=M_{\textrm{Z},{\star}}/M_{\star}\simeq10^{-4}$
while
$M_{\textrm{g}}\gt0$
.
For comparison, the behaviour in the
$M_{\textrm{h}} = 10^{10}\,\text{M}_{\odot}$
mass bin (lower panels) reveals that the mode of star formation – instantaneous or bursty – has little-to-no effect on the growth of
$M_{\textrm{g}}$
,
$M_{\textrm{Z},\textrm{g}}$
,
$M_{\star}$
or
$M_{\textrm{Z},{\star}}$
. Given this lack of sensitivity to star formation mode at higher masses, we focus on the
$M_{\textrm{h}} = 10^7\,\text{M}_{\odot}$
mass bin in the remainder of this subsection.
Having established the fiducial behaviour of the instantaneous and delayed feedback scenarios, for the rest of this work, we limit our focus to only the delayed feedback scenario (but see Section 3.3).
Sensitivity to IGM metallicity,
${\boldsymbol{Z}}_{\textbf{{IGM}}}$
: Figure 2 shows how our model predictions depend on what we assume for the metallicity of IGM,
$Z_{\text{IGM}}$
. This sets the metallicity of the cosmologically accreted gas (cf. the first term in Equation 4). The upper (lower) panel shows the case for
$Z_{\text{IGM}}$
=
$10^{-5}\,\text{Z}_{\odot}$
(
$10^{-1}\,\text{Z}_{\odot}$
); recall that the fiducial value is
$10^{-3}\,\text{Z}_{\odot}$
(grey bands and curves correspond to
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
for this fiducial case). Note that we do not expect the value of
$Z_{\text{IGM}}$
to affect the values of
$M_{\textrm{g}}$
and
$M_{\star}$
in any way and hence we do not show the interdecile ranges for these quantities.
Influence of the IGM metallicity,
${\textbf{{Z}}}_{\text{IGM}}$
: We show the evolution of
$M_{\textrm{g}}$
,
$M_{\star}$
,
$M_{\textrm{Z},\textrm{g}}$
, and
$M_{\textrm{Z},{\star}}$
(solid, dot-dashed, dashed, dotted curves, respectively) with cosmic time/redshift as the metallicity of the accreted gas
$(Z_{\mathrm{IGM}})$
is varied. The upper and lower panels, respectively, corresponds to
$Z_{\text{IGM}}=10^{-5}\,\text{Z}_{\odot}$
and
$10^{-1}\,\text{Z}_{\odot}$
. Grey bands and curves correspond to
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
for the fiducial
$Z_{\text{IGM}}$
.
For
$Z_{\text{IGM}}=10^{-5}\,\text{Z}_{\odot}$
(upper panels), the differences with respect to the fiducial predictions are negligible except at early cosmic times –
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
are marginally lower than in the fiducial case by less than 0.1 dex at
$z\gtrsim 11$
. For
$Z_{\text{IGM}}=10^{-1}\,\text{Z}_{\odot}$
(lower panels), differences with respect to fiducial predictions are more readily apparent – both
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
track the fiducial case but are offset by approximately 1 dex.
Influence of the heavy element yield,
${{{Y}}}_{{\textrm{Z}}}$
: We show the evolution of
$M_{\textrm{g}}$
,
$M_{\star}$
,
$M_{\textrm{Z},\textrm{g}}$
, and
$M_{\textrm{Z},{\star}}$
(solid, dot-dashed, dashed, dotted curves, respectively) with cosmic time/redshift as heavy elements’ yield
$(Y_{\textrm{Z}})$
is varied. Upper and lower panels corresponds to
$Y_{\textrm{Z}}=0.006$
and
$Y_{\textrm{Z}}=0.6$
. Grey bands and curves correspond to
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
for the fiducial
$Y_{\textrm{Z}}$
.
Influence of metal-dependent feedback: We show the impact of metal-dependent feedback (Equation 8) on the evolution of
$M_{\textrm{g}}$
,
$M_{\star}$
,
$M_{\textrm{Z},\textrm{g}}$
, and
$M_{\textrm{Z},{\star}}$
(solid, dot-dashed, dashed, dotted curves, respectively) for lower and higher halo masses at
$z=5$
. The upper and lower panels correspond to
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
and
$10^{10}\,\text{M}_{\odot}$
, respectively. As before, the grey curves correspond to the fiducial model.
Sensitivity to heavy element yield,
${\textbf{{Y}}}_{\textbf{{Z}}}$
: Figure 2 shows how our model predictions are affected by what we assume for heavy element yield,
$Y_{\textrm{Z}}$
. This tracks the production of metals by (primarily massive) stars, which proceed to enrich the gas at the end of their main sequence life and are released into the galaxy’s gaseous reservoir. The upper (lower) panel shows the case for
$Y_{\textrm{Z}}=0.006\,(0.6)$
; recall that the fiducial value is 0.06. As in Figure 2, grey bands and curves correspond to
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
for this fiducial value of
$Y_{\textrm{Z}}$
. As for variations in
$Z_{\text{IGM}}$
, we do not expect the value of
$Y_{\textrm{Z}}$
to affect the values of
$M_{\textrm{g}}$
and
$M_{\star}$
.
For
$Y_{\textrm{Z}}=0.006$
(upper panel) and
$Y_{\textrm{Z}}=0.6$
(lower panel), the differences with respect to the fiducial predictions are straightforward – a decreased (increased) heavy element yield results in decreased (increased) values of
$M_{\textrm{Z},\textrm{g}}$
and
$M_{\textrm{Z},{\star}}$
. The evolution of these quantities follows that in the fiducial case, with similar shapes but offset by
$\sim$
1 dex below (above) the fiducial curves for
$Y_{\textrm{Z}}=0.006$
(0.6).
Metal-dependent feedback: As noted earlier, there are good physical reasons to expect that the efficiency of feedback should be metal-dependent – for example, because higher mass stars collapse directly into black holes (e.g. O’Connor & Ott Reference O’Connor and Ott2011; Sukhbold et al. Reference Sukhbold, Ertl, Woosley, Brown and Janka2016; Jecmen & Oey Reference Jecmen and Oey2023) or because the absence of metals in their outer envelopes reduces the mass and momentum flux of winds (cf. Lamers & Cassinelli 1999). Figure 4 shows the impact of our assumed form for metal-dependent feedback (cf. Section 2.2.3) on our model predictions for lower and higher halo masses –
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
(upper panel) and
$10^{10}\,\text{M}_{\odot}$
(lower panel). Equation (8) for
$f(Z_{\star})$
implies that there is a decrease in total momentum at sub-solar metallicities, and therefore a decrease in feedback efficiency at lower metallicities.
The effect of the reduced feedback efficiency is readily evident for the case of
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
: there is stronger growth in the gaseous and stellar components when the feedback efficiency is metal-dependent:
$M_{\star}$
is offset by
$\sim\!0.5$
dex,
$M_{\textrm{Z},{\star}}$
by
$\sim\!1$
dex.
The differences between
$M_{\textrm{g}}$
and
$M_{\textrm{Z},\textrm{g}}$
, relative to the fiducial case, are more marked and qualitative. The weaker feedback at low metallicities and early times acts to erase the oscillations that are so apparent in the fiducial case.
$M_{\textrm{g}}$
has an amplitude similar to that in the fiducial case, but the median system can retain the gas mass for
$\sim$
$0.1$
Gyr longer. The same behaviour is evident in the median value of
$M_{\textrm{Z},\textrm{g}}$
, while the large fluctuations in the 10th-to-90th variation evident in the fiducial case are eliminated.
This contrasts with what we observe at the higher halo mass of
$10^{10}\,\text{M}_{\odot}$
, in the lower panel. At these masses, it is evident that the impact of metal-dependent feedback is negligible. The feedback mechanism plays a secondary role compared to the dominant influence of the deep gravitational potential of these haloes.
3.2. Mass-metallicity relations
In this subsection, we turn our attention to our model predictions for the relationship between stellar mass and the median gas-phase (
$Z_g=M_{\textrm{Z},\textrm{g}}/M_{\textrm{g}}$
) and stellar metallicities (
$Z_{\ast}=M_{\textrm{Z},{\star}}/ M_{\star}$
) at different epochs – Figures 5 and 6, respectively. As in the previous subsection, we focus on results for the delayed feedback scenario.
Gas-phase metallicity: In Figure 5, we show how gas-phase metallicity,
$Z_{\textrm{g}}$
, varies with stellar mass,
$M_{\star}$
, at redshifts
$z=5$
, 7, and 10 (blue, green and red curves, respectively), for our fiducial feedback model (left panel) and metal-dependent feedback (right panel). The curves represent the median behaviour for a given
$M_{\star}$
at that epoch; the associated coloured bands indicate the
$10^{\textrm{th}}$
to
$90^{\textrm{th}}$
percentile variation. The symbols correspond to observational data from Arellano-Córdova et al. (Reference Arellano-Córdova2022), Nakajima et al. (Reference Nakajima2023), Trump et al. (Reference Trump2023), Chemerynska et al. (Reference Chemerynska2024), and Curti et al. (Reference Curti2024).
The fiducial model predicts that the metallicity at higher stellar mass (
$M_{\star}$
) is relatively flat at a given epoch, with
$Z_{\textrm{g}} \simeq 0.03\,\text{Z}_{\odot}$
. However, at lower masses, the metallicity declines sharply, with a corresponding increase in the scatter in
$Z_{\textrm{g}}$
. The mass scale at which this decline occurs increases with decreasing redshift. At
$z=7$
and
$z=10$
, there is a monotonic decline in
$Z_{\textrm{g}}$
. However, at
$z=5$
, we observe a sharp upturn followed by a precipitous drop. This behaviour at
$z=5$
is an artifact of the onset of UV suppression in
$z_{\text{rei}}=7$
, which introduces a large scatter in the values of
$M_{\textrm{g}}$
and
$M_{\textrm{Z},\textrm{g}}$
. This scatter is evident in Figure 1, particularly in the
$10{\text{th}}$
to
$90{\text{th}}$
percentile variation, and affects
$Z_{\textrm{g}}$
due to the small value ratios. Some of this scatter arises from stochastic halo growth histories in the merger trees, which, together with UV suppression, amplify the apparent upturn and subsequent drop.
In contrast, the effect of metal-dependent feedback is to suppress the decrease in
$Z_{\textrm{g}}$
at low
$M_{\star}$
and
$M_{\textrm{h}}$
at
$z=7$
and
$z=10$
. The relation remains relatively flat, with
$Z_{\textrm{g}} \simeq 0.03$
–
$0.04\,\text{Z}_{\odot}$
across the entire mass range. A decline in
$Z_{\textrm{g}}$
becomes evident only at
$z=5$
, showing the same sharp upturn and precipitous drop seen in the instantaneous case, but occurring at a stellar mass roughly a factor of
$\sim$
10 smaller.
The decline in gas-phase metallicity (
$Z_{\textrm{g}}$
) with decreasing mass at
$z = 10$
and
$z = 7$
reflects the increasing efficiency of feedback-driven metal ejection in shallow potential wells. By
$z = 5$
, this trend reverses: low-mass halos exhibit higher
$Z_{\textrm{g}}$
because reionisation-driven quenching has already expelled most of their gas, and metal-enriched ejecta from the most recently formed stars enrich the small remaining reservoirs. This process produces the observed upturn in our model. At even lower masses, stellar feedback combined with UV suppression of cosmological gas accretion is sufficient to remove gas entirely from these systems.
These predictions imply that multi-band photometry combined with metallicity measurements can discriminate between feedback scenarios. For example, delayed feedback generates larger scatter in the
$Z_{\textrm{g}}$
–
$M_{\star}$
relation (
$\Delta\log Z_{\textrm{g}} \sim 0.5$
–1 dex) and sharper declines at characteristic masses that evolve with redshift.
Stellar metallicity: Figure 6 mirrors Figure 5; here we show how stellar metallicity,
$Z_{\star}$
, varies with stellar mass,
$M_{\star}$
at different redshifts assuming our fiducial and metal-dependent feedback models (left and right panels respectively).
We see that
$Z_{\star}$
decreases with decreasing
$M_{\star}$
for all redshifts. There is strong overlap between the relations at different redshifts, as we might expect;
$M_{\star}$
and
$M_{\textrm{Z},{\star}}$
will either remain constant or grow in time, which contrasts with
$M_{\textrm{g}}$
and
$M_{\textrm{Z},\textrm{g}}$
, which can either increase or decrease in response to cosmological gas accretion, enrichment episodes via supernovae, and expulsion via stellar-driven winds. Compared to the behaviour we see in
$Z_{\textrm{g}}$
with
$M_{\ast}$
, there is a sharper decline in
$Z_{\ast}$
with decreasing
$M_{\ast}$
and an increase in the 10th–90th percentile variation.
Stellar mass versus gas-phase metallicity relations between
$z=5$
and 10: Here we plot gas phase metallicity,
$Z_{\textrm{g}}=M_{\textrm{Z},\textrm{g}}/M_{\textrm{g}}$
, as a function of stellar mass,
$M_{\star}$
, in units of
$\text{M}_{\odot}$
at
$z=5$
, 7, and 10 (blue, green and red curves respectively). The left panel shows the behaviour in our fiducial feedback model; the right shows the impact of our assumed metal-dependent feedback model. Filled symbols correspond to observational data from Arellano-Córdova et al. (Reference Arellano-Córdova2022), Nakajima et al. (Reference Nakajima2023), Trump et al. (Reference Trump2023), Chemerynska et al. (Reference Chemerynska2024), and Curti et al. (Reference Curti2024).
Stellar mass versus stellar metallicity relations between
$z=5$
and 10: We show the same information as in Figure 5, but for stellar metallicity,
$Z_{\ast}=M_{\textrm{Z},{\star}}/M_{\star}$
, rather than gas-phase metallicity.
3.3. Evolution of dust,
$M_{\textrm{d}}$
In this final subsection, we examine our model predictions for the dust mass,
$M_{\textrm{d}}$
. As usual, we focus on two representative systems at
$z=5$
: low-mass halos with
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
and high-mass halos with
$M_{\textrm{h}}=10^{10}\,\text{M}_{\odot}$
, and we consider the behaviour expected for both instantaneous and delayed feedback scenarios. Figures 7 and 8 show the evolution of gas, stellar, and dust masses (solid, dot-dashed, and dashed curves, respectively) as a function of cosmic time for these two halo masses. As in our other plots, we follow both the median trend and the
$10^{\textrm{th}}$
–
$90^{\textrm{th}}$
percentile range of
$M_{\textrm{d}}$
.
Across both mass scales,
$M_{\textrm{d}}$
broadly traces the behaviour of
$M_{\textrm{g}}$
, offset in normalisation by a factor of
$\sim\!1/100$
, with a similar level of scatter. This is expected: our model links dust growth to the SFR,
$\dot{M}_{\ast}$
, which depends directly on
$M_{\textrm{g}}$
. The evolution also mirrors that of the gas-phase metal mass
$M_{\textrm{Z},\textrm{g}}$
(see Figure 1), reflecting the fact that, to the first order,
$M_{\textrm{d}}$
can be obtained from
$M_{\textrm{Z},\textrm{g}}$
via a constant metal-to-dust ratio. Consequently, the oscillatory behaviour of
$M_{\textrm{Z},\textrm{g}}$
is echoed in
$M_{\textrm{d}}$
. In our model, deviations from this constant ratio occur whenever dust destruction channels (Equations 11–15) operate more efficiently than metal loss, particularly in high-mass halos where supernovae shocks dominate. The scatter in
$M_{\textrm{d}}$
for a given SFR – which at fixed halo mass
$M_{\textrm{h}}$
will show a scatter as a result of variations in assembly history – implies that dust-based SFR estimates could be significantly biased if this intrinsic variation is not taken into account.
For low-mass
$M_{\textrm{h}}=10^7\,\text{M}_{\odot}$
halos,
$M_{\textrm{g}}$
almost never exceeds
$10^5\,\text{M}_{\odot}$
, so the primary pathway for dust removal is through stellar winds; these systems never reach the Sedov–Taylor phase of supernova evolution. In the instantaneous feedback scenario (upper panel of Figure 7),
$M_{\textrm{d}}$
follows the gradual decline of
$M_{\textrm{g}}$
once the UVB is turned on at
$z_{\mathrm{rei}}$
. In contrast, in the delayed feedback scenario (lower panel),
$M_{\textrm{d}}$
declines to zero along with
$M_{\textrm{g}}$
, as the combined effects of the UVB and feedback rapidly reduce the gas reservoirs in these galaxies.
The high-mass
$M_{\textrm{h}}=10^{10}\,\text{M}_{\odot}$
halos present a different picture. Their deeper gravitational potential wells allow them to retain more gas, enabling supernova shocks – in addition to stellar winds – to destroy dust. Once
$M_{\textrm{g}}$
is high enough for the Sedov–Taylor phase to develop, dust destruction becomes the dominant removal process, producing sharp declines in
$M_{\textrm{d}}$
. Only when
$M_{\textrm{g}} \gtrsim 10^7\,\text{M}_{\odot}$
does dust production begin to balance destruction. As a result, our model predicts that high-mass systems can host ongoing star formation while exhibiting a wide diversity in dust content, from almost dust-free to dust-rich states. This diversity is qualitatively consistent with the oscillatory behaviour seen in more sophisticated simulations (e.g. Choban et al. Reference Choban, Salim, Kereš, Hayward and Sandstrom2025) and could reconcile observations of high-z galaxies with a wide range of dust masses (e.g. Rodighiero et al. Reference Rodighiero2023; Barrufet et al. Reference Barrufet2023)
Evolution of dust mass
$M_{\textrm{d}}$
for low-mass halos: The growth of
$M_{\textrm{g}}$
(solid),
$M_{\star}$
(dash-dotted),
$M_{\textrm{d}}$
(dashed) in halos of mass bin (at
$z=5$
)
$M_{\textrm{h}}=10^{7}\,\text{M}_{\odot}$
. The upper(lower) panel is for the instantaneous(delayed) feedback scenario. For clarity we only show the median behaviour for
$M_{\textrm{g}}$
and
$M_{\star}$
while the shaded region represents the
$10\mathrm{th}$
–
$90\mathrm{th}$
variation in
$M_{\textrm{d}}$
.
Evolution of dust mass
$M_{\textrm{d}}$
for high-mass halos: Similar to Figure 7 we show the
$M_{\textrm{g}}$
(solid),
$M_{\star}$
(dash-dotted),
$M_{\textrm{d}}$
(dashed) but for halos of mass bin (at
$z=5$
)
$M_{\textrm{h}}=10^{10}\,\text{M}_{\odot}$
.
4. Discussion
We have demonstrated that bursty star formation fundamentally alters the chemical evolution of early galaxies compared to continuous star formation scenarios. The delayed feedback inherent to bursty star formation drives oscillatory behaviour in both gas-phase and stellar metallicities, while reionisation introduces additional complexity by suppressing gas accretion. Together, these processes generate large scatter in gas mass (
$M_{\textrm{g}}$
) and gas-phase metallicity (
$M_{\textrm{Z},\textrm{g}}$
), producing distinct signatures in the mass–metallicity relation that depend on the feedback prescription.
Our predictions align closely with recent JWST observations, particularly for galaxies with
$M_{\star} \sim 10^6 {-} 10^8$
M
$_{\odot}$
. Figure 5 shows that delayed feedback combined with UV suppression reproduces the low normalisation (
$M_{\textrm{Z},\textrm{g}} \sim 0.02{-}0.05\,Z_{\odot}$
) and large intrinsic scatter (
$\Delta \log M_{\textrm{Z},\textrm{g}} \sim 0.5{-}1$
dex) observed at
$z \simeq 5{-}10$
. The steep decline in
$M_{\textrm{Z},\textrm{g}}$
at low masses and its evolution with redshift match the trends reported by Nakajima et al. (Reference Nakajima2023), Curti et al. (Reference Curti2024), and Chemerynska et al. (Reference Chemerynska2024), while the upturn and subsequent drop at
$z \simeq 5$
reflect reionisation-driven gas depletion inflating
$M_{\textrm{Z},\textrm{g}}$
in residual reservoirs before complete gas removal. Metal-dependent feedback flattens the relation at
$z \simeq 7\!-\!10$
, yielding
$M_{\textrm{Z},\textrm{g}} \sim 0.03\!-\!0.04\,Z_{\odot}$
across the mass range, consistent with the lower end of the observed normalisation.
Observational metallicities from Arellano-Córdova et al. (Reference Arellano-Córdova2022) and Sarkar et al. (Reference Sarkar2025) agree with the plateau and sharp decline predicted by our delayed feedback models at
$z \simeq 6\!-\!10$
. Trump et al. (Reference Trump2023) and Li et al. (Reference Li2025) report higher values – 12+
$\log(\mathrm{O/H}) \simeq 7.2\!-\!8.0$
(
$\sim 0.03\!-\!0.18\,Z_{\odot}$
) – that more closely resemble our instantaneous feedback + UV background results, highlighting the observational split between feedback scenarios. The lower normalisation of the mass–metallicity relation at
$z \simeq 7{-}10$
noted by Heintz et al. (Reference Heintz2022) and Ucci et al. (Reference Ucci2022) is also reproduced, suggestive of pristine inflows and strong outflows in low-mass halos. Furthermore, the decoupling of gas-phase and stellar metallicities in our models mirrors episodic enrichment histories inferred by Rey et al. (Reference Rey2025) and elevated N/O ratios reported by Ji et al. (Reference Ji2025).
Theoretical and simulation studies reinforce these trends. The DELPHI23 semi-analytic model (Mauerhofer et al. Reference Mauerhofer2025) predicts a sharp metallicity drop at low masses, similar to our delayed + UV results at
$z \simeq 5$
. Cosmological hydrodynamical simulations such as FIRE and FIRE-2 (Ma et al. Reference Ma2015; Hopkins et al. Reference Hopkins2018) exhibit mass–metallicity relations whose shapes and scatter closely match our predictions for
$M_{\star} \sim 10^6{-}10^8\,\text{M}_{\odot}$
, including steep declines in
$M_{\textrm{Z},\textrm{g}}$
in burst-dominated regimes. In FIRE, the stellar metallicity–stellar mass relation also resembles our delayed + UV case, with offsets of only
$\sim 0.5{-}2$
dex at comparable redshifts.
Our dust evolution predictions capture key features seen in recent high-z surveys. The delayed growth and large scatter in dust mass (
$M_{\textrm{d}}$
) in our models parallel the diversity observed by Rodighiero et al. (Reference Rodighiero2023) and Barrufet et al. (Reference Barrufet2023), and match the quenching or decline in dust content seen in FIRE-2 simulations (Choban et al. Reference Choban, Salim, Kereš, Hayward and Sandstrom2025). In contrast, the CROC simulations Esmerian & Gnedin (Reference Esmerian and Gnedin2022) predict steady dust growth in massive halos, whereas our delayed feedback scenario produces a pronounced quench over similar timescales - likely reflecting our assumption that
$M_{\mathrm{ISM}} = M_{\textrm{g}}$
, which may overestimate dust destruction in gas-rich halos.
Taken together, these comparisons show that burst-driven, delayed feedback, when combined with reionisation suppression, can naturally explain several observed high-z trends – oscillatory enrichment, large scatter in
$M_{\textrm{Z},\textrm{g}}$
, decoupled stellar and gas-phase metallicities, and stochastic dust evolution. At the same time, the split between observational cases favouring delayed versus instantaneous feedback underscores the need for joint constraints using metallicity, dust content, and star formation variability to fully pin down feedback timescales in early galaxies.
5. Conclusions
We have investigated how bursty star formation during cosmic reionisation affects the metallicity and dust content of low-mass galaxies using our Ashvini model. Bursty star formation can drive strong outflows, which suppress accretion of pristine gas from the intergalactic medium. This can offset metallicity in these galaxies, and the episodic nature of enrichment can decouple gas-phase and stellar metallicities, which will influence how we interpret the observable properties of the high-redshift galaxy population. Our results confirm these expectations:
-
• Bursty star formation, driven by delayed feedback, induces oscillations in gas and stellar metallicities (
$Z_{\textrm{g}}$
,
$Z_{\star}$
), particularly in lower-mass halos, as illustrated by our
$M_{\textrm{h}} \sim 10^7\,\text{M}_{\odot}$
example (Figure 1). This leads to a decoupling of
$Z_{\textrm{g}}$
and
$Z_{\star}$
, consistent with the episodic accretion and enrichment histories inferred empirically by Chemerynska et al. (Reference Chemerynska2024) using JWST data and theoretically by Rey et al. (Reference Rey2025) using the EDGE simulations (Engineering Dwarfs at Galaxy Formation’s Edge; cf. Agertz et al. Reference Agertz2020) of dwarf galaxies. -
• Reionisation introduces significant scatter in
$Z_{\textrm{g}}$
, producing a sharp upturn followed by a decline at
$z \lesssim 7$
(Figure 1, lower panel). This feature aligns with the metallicity plateau and subsequent drop observed in JWST studies (e.g. Arellano-Córdova et al. Reference Arellano-Córdova2022; Sarkar et al. Reference Sarkar2025). -
• Metal-dependent feedback flattens the
$Z_{\textrm{g}}$
–
$M_{\star}$
relation to
$Z_{\textrm{g}} \approx 0.03$
–
$0.04\,\text{Z}_{\odot}$
at
$z = 7$
–10 (Figure 5), in agreement with the lower normalisation of the mass–metallicity relation reported by Heintz et al. (Reference Heintz2022) and Ucci et al. (Reference Ucci2022). -
• Dust evolution broadly tracks
$Z_{\textrm{g}}$
but is delayed by feedback cycles. The resulting diversity in dust content, especially in high-mass halos, is consistent with ALMA and JWST observations of dust-rich and dust-poor galaxies at
$z \gtrsim 6$
Rodighiero et al. (Reference Rodighiero2023) and Barrufet et al. (Reference Barrufet2023).
The results confirm that inferences of dust content or stellar age based on galaxy colour may be misleading in bursty, metal-poor systems. Because dust production follows gas-phase metallicity and is delayed by feedback, both the timing and amount of dust are sensitive to burst histories. Bursty star formation introduces time lags and variability in enrichment. The metallicity of gas inflow properties and the metallicity-dependence of supernovae can critically shape metallicity and dust evolution. These effects are especially strong in low-mass galaxies, where shallow potentials amplify burst effects. As such, observations of galaxies in the early Universe must account for these burst-driven dynamics when interpreting metallicity, dust, and star formation rates.
Acknowledgements
The authors thank the anonymous referee for their careful and thoughtful report. SB acknowledges the Dr Albert Shimmins Fund for a Postgraduate Writing-Up Award, and support from grant PID2022-138855NB-C32 funded by MICIU/AEI/10.13039/501100011033 and ERDF/EU, and project PPIT2024-31833, cofunded by EU-Ministerio de Hacienda y Función Pública-Fondos Europeos-Junta de Andaluca-Consejera de Universidad, Investigación e Innovación. CP acknowledges the support of the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. This research uses python (van Rossum Reference van Rossum1995) libraries, including numpy (Harris et al. Reference Harris2020), matplotlib (Hunter Reference Hunter2007), and scipy (Virtanen et al. Reference Virtanen2020).
Code and data availability
The Ashvini model – along with the dark matter halo merger trees used in this work – is made publicly available at https://github.com/Anand-JM97/Ashvini.
Appendix A. Calculation of
$\varepsilon_{\text{uv}}$
In Reference Menon and PowerMP24, we used the parameter
$\epsilon_{\text{in}}\equiv\,\varepsilon_{\text{UV}}$
to regulate cosmological gas accretion onto halos in the presence of a UVB, arguing that suppression of accretion will be dominated by the UVB at redshifts
$z\gtrsim 5$
. For completeness we summarise the arguments of Section 2.4.1 in Reference Menon and PowerMP24, which builds on Gnedin (Reference Gnedin2000), Okamoto, Gao, & Theuns (Reference Okamoto, Gao and Theuns2008), and Kravtsov & Manwadkar (Reference Kravtsov and Manwadkar2022), to motivate the functional form of
$\varepsilon_{\text{UV}}$
.
The baryon fraction within a halo is given by
\begin{equation} f_{\textrm{b}}(M_{\textrm{h}},z)=\frac{\Omega_{\textrm{b}}}{\Omega_{\textrm{m}}}s(\mu_{\textrm{c}},\omega)\end{equation}
where
$\Omega_{\textrm{b}}/\Omega_{\textrm{m}}$
is the cosmic baryon fraction. s(x, y) is defined as,
\begin{equation} s(x,y)= \big[1+(2^{\frac{y}{3}}-1)x^{-\textrm{y}}\big]^{-\frac{3}{y}},\end{equation}
where
$\mu_{\textrm{c}}=m_{\textrm{h}}/M_{\textrm{c}}(z)$
.
$M_{\textrm{c}}(z)$
is a characteristic mass below which
$s(\mu_{\textrm{c}},\omega) \rightarrow 0$
,
\begin{align} M_{\textrm{c}}=1.69\times 10^{10}\frac{\text{exp}{(\!-0.63z)}}{1+\text{exp}([z/\beta]^{\zeta})}M_{\odot}, \end{align}
where
$\zeta=15$
and
$\beta$
is given by,
\begin{equation} \beta=z_{\text{rei}}\left[\text{ln}(1.82\times10^3\text{exp}{(\!-0.63z_{\text{rei}})}-1)\right]^{-1/\zeta}\! ;\end{equation}
here
$z_{\text{rei}}$
is the redshift of reionisation.
We write,
\begin{align}\epsilon_{\text{in}}=\text{max}\bigg(0,s(\mu_{\textrm{c}},\omega)\bigg[(1+X)- 2\epsilon(z,\zeta)\frac{M_{\textrm{h}}}{\dot{M}_{\textrm{h}}}X(1+z)H(z)\bigg]\bigg), \end{align}
where
\begin{align*}\epsilon(z,\zeta)&=\frac{0.63}{1+e^{(z/\beta)^\zeta}}+\frac{\zeta z^{\zeta-1}}{\beta^\zeta}\frac{e^{(z/\beta)^\zeta}}{(1+e^{(z/\beta)^\zeta})^2},\\[9pt]X&=\frac{3c_{\omega} M_{\omega}}{1+c_{\omega} M_{\omega}},\quad c_{\omega}=2^{\omega/3}-1,\\[9pt]M_{\omega}&=\left(\frac{M_{\textrm{c}}(z)}{M_{\textrm{h}}}\right)^\omega,\quad\omega=2.\end{align*}
The net effect is to suppress gas accretion at halo masses
$M_{\textrm{h}}\leq\!M_{\textrm{c}}(z)$
at
$z \lt z_{\text{rei}}$
, with
$M_{\textrm{c}}(z)$
increasing with decreasing z (i.e. accretion is suppressed onto progressively higher mass halos). For redshifts
$z \gt z_{\text{rei}}$
, there is no suppression of gas accretion.
Open access
