The Cosmic History of Long Gamma-Ray Bursts

The distinctive features of long-duration gamma-ray bursts (LGRBs)—such as their detection up to high redshifts (z = 8.2–9.2; Salvaterra et al. 2009; Tanvir et al. 2009; Cucchiara et al. 2011), their association with faint blue star-forming galaxies, their positional coincidence with the brightest ultraviolet regions within their hosts (Fruchter et al. 2006), and their association with core-collapse supernovae (Hjorth & Bloom 2012)—point to a massive-star progenitor, thus establishing a direct link with the formation and relatively rapid evolution of massive stars across cosmic times. The study of increasingly larger samples of hosts has revealed the preference of LGRBs to form in subsolar-metallicity environments (Vergani et al. 2015; Perley et al. 2016; Palmerio et al. 2019). This is in agreement with the progenitor model (MacFadyen & Woosley 1999) requiring a low metallicity to maintain a high specific angular momentum at the time of the core collapse, in order to produce the GRB jet (Yoon et al. 2006). Therefore, the physical conditions for the production of a relativistic jet determine the intrinsic properties of the LGRB population, such as their characteristic luminosities and cosmic rates.

GRB population studies aim at reconstructing the population luminosity and formation rate functions. The free parameters of these functions are often constrained by reproducing sets of rest-frame and observer-frame properties (Daigne et al. 2006; Wanderman & Piran 2010; Robertson & Ellis 2012) of the detected GRB samples and properly accounting for sample incompleteness (Salvaterra et al. 2012; Ghirlanda et al. 2015; Palmerio & Daigne 2021).

In this work, we build a GRB population model (Section 2) and derive the cosmic formation rate of LGRBs and their luminosity function, together with its evolution with redshift. To this end, we exploit, for the first time, the largest set of observational constraints collected by different GRB detectors over the last four decades (Section 3). We account for the presence of relativistic jets and the orientation-dependence of the observed GRB properties. We derive the evolution of the ratio of the GRB formation rate to the cosmic star formation rate (CSFR), and test its consistency with the host metallicity bias (Section 4). We discuss the implications of our population model in Section 5. Throughout the paper, we assume a standard flat cosmology, with Ω_M = 0.3 and h = 0.7.

We build a synthetic LGRB population model, based on two distribution functions. The GRB formation rate per unit of comoving volume ρ(z) is modeled with the parametric function (Cole et al. 2001; Madau & Dickinson 2014)

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}\displaystyle \frac{{\left(1+z\right)}^{{p}_{z,1}}}{1+{\left(\tfrac{1+z}{{p}_{z,2}}\right)}^{{p}_{z,3}}},\end{eqnarray} \tag{ 1 }$$

where ρ₀, in Gpc⁻³ yr⁻¹, is the LGRB local density rate. The choice of this functional form is motivated by its relatively small number of free parameters and, for ease of comparison, by its systematic use for fitting the CSFR (Madau & Dickinson 2014; Madau & Fragos 2017; Hopkins & Beacom 2006). Flat priors are considered for the free parameters (p_z,1, p_z,2, and p_z,3) that can vary over wide ranges, also including the values corresponding to the CSFR. This ensures that we can recover, if present, any significant deviation of the GRB formation rate from the CSFR.

The luminosity function of LGRBs is often described in the literature by a broken power law (e.g., Salvaterra et al. 2009; Salvaterra et al. 2012):

$$\begin{eqnarray}{\rm{\Phi }}({L}_{\mathrm{iso}},z)\propto \left\{\begin{array}{ll}{L}_{\mathrm{iso}}^{-{a}_{1}} & \mathrm{if}\,\,{L}_{\mathrm{iso}}\leqslant {L}_{\mathrm{iso},{\rm{b}}}={L}_{\mathrm{iso},0}{\left(1+z\right)}^{\delta }\\ {L}_{\mathrm{iso}}^{-{a}_{2}} & \mathrm{if}\,\,{L}_{\mathrm{iso}}\gt {L}_{\mathrm{iso},{\rm{b}}}={L}_{\mathrm{iso},0}{\left(1+z\right)}^{\delta }\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where a₁ and a₂, respectively, are the indices of the power laws below and above the break luminosity L_iso,b, which is left free to evolve with redshift ∝(1 + z)^δ . Therefore, L_iso,0 represents the break of the luminosity function at z = 0. Φ(L_iso, z) is normalized to its integral and it is defined for L_iso ≥ 10⁴⁷ erg s⁻¹. An alternative working hypothesis sometimes considered in the literature is the Schechter (Schechter 1976) luminosity function. We discuss in Appendix A.3 the results obtained under this different assumption. In summary, our model allows for the evolution of the GRB population both in terms of rate density (Equation (1)) and luminosity (Equation (2)).

LGRBs follow the E_p − L_iso (Yonetoku et al. 2004) correlation, where E_p is the peak energy of the ν F_ν spectrum and the E_p − E_iso correlation (Amati & Frontera 2002) links the peak energy with the isotropic energy E_iso. It has been proved (Nava et al. 2012; Dai et al. 2021) that these correlations do not evolve with redshift. Following Ghirlanda et al. (2015), in our simulation we sample Equation (1), extracting a redshift, and we independently assign an E_p value sampling a broken power law. The values of L_iso and E_iso corresponding to E_p are derived from the E_p − L_iso and E_p − E_iso correlations, respectively, accounting for their Gaussian scatters (Nava et al. 2012). For ease of comparison with the literature, once we have derived the posterior parameter distributions of the model, we map, through the E_p − L_iso correlation, the peak energy distribution function back to the luminosity function, and we provide the resulting parameter values in Table 1. GRB rest-frame durations are estimated as f · E_iso/L_iso, where f is a scaling factor sampled from a free Gaussian distribution ${\mathfrak{G}}[f| \mu ,\sigma ]$.

Table 1. Parameters of the GRB Population Model

ρ(z)	${p}_{z,1}={3.33}_{-0.33}^{+0.33}$	${p}_{z,2}={3.42}_{-0.28}^{+0.28}$	${p}_{z,3}={6.21}_{-0.32}^{+0.38}$	${\rho }_{0}={79}_{-33}^{+57}$ Gpc−3 yr−1	⋯
Φ(Liso, z)	${a}_{1}\,=\,{0.97}_{-0.04}^{+0.05}$	${a}_{2}\,=\,{2.21}_{-0.18}^{+0.13}$	$\mathrm{log}{L}_{\mathrm{iso},0}={52.02}_{-0.19}^{+0.22}$	$\delta ={0.64}_{-0.26}^{+0.32}$	⋯
${\mathfrak{G}}[\mathrm{log}f| \mu ,\sigma ]$	$\mu ={3.00}_{-0.26}^{+0.29}$	$\sigma ={0.23}_{-0.04}^{+0.04}$	⋯	⋯	⋯
${\mathfrak{G}}[\alpha | \mu ,\sigma ]$	μ = −1.0	σ = 0.1	⋯	⋯	⋯
${\mathfrak{G}}[\beta | \mu ,\sigma ]$	μ = −2.3	σ = 0.1	⋯	⋯	⋯
	K	S	Σ	Y0	X0
$\mathrm{log}$ Eiso- $\mathrm{log}$ θjet	${9.65}_{-0.01}^{+0.01}$	$-{0.203}_{-0.002}^{+0.002}$	0.1	1	1
$\mathrm{log}$ Ep- $\mathrm{log}$ Liso	0.27	0.56	0.18	406.7	3.38 × 1052
$\mathrm{log}$ Ep- $\mathrm{log}$ Eiso	−0.12	0.58	0.20	406.7	9.10 × 1052
$\mathrm{log}$Γ0- $\mathrm{log}$ Eiso	2.1	0.3	0.1	1	1052

Note. The parameters constrained through the MCMC (Section 2) are reported with the associated 1σ uncertainty. Their posterior distributions are shown in Figure 2. The parameters without errors are held constant. The four correlations listed are expressed in the form $\mathrm{log}\left(\displaystyle \frac{Y}{{Y}_{0}}\right)=K+S\cdot \mathrm{log}\left(\tfrac{X}{{X}_{0}}\right)$, with the parameter Σ representing their scatter along the independent coordinate. E_iso, L_iso, E_p, and θ_jet have units of erg, erg s⁻¹, keV, and rad, respectively.

Download table as:  ASCIITypeset image

It is known that GRBs are highly collimated (e.g., Ghirlanda et al. 2007). Therefore, we compute their intrinsic rate ρ(z) (Equation (1)) by accounting for the elusive bulk of the population consisting of events with jets oriented to the side of our line of sight. We assume uniform conical jets distributed with random orientations θ_v with respect to the line of sight. The jet opening angle values θ_jet appear to be correlated with E_iso (Ghirlanda et al. 2007), such that we assume a free parametric correlation θ_jet–E_iso. Even off-axis jets (θ_v > θ_jet) may be detected, depending on their orientation and bulk Lorentz factor Γ, characterizing the prompt emission phase. We account, for the first time within a population study of LGRBs, for the relativistic effects on the observed flux and fluence, as well as the peak energy induced by the jet orientation and the beaming through the Doppler factor $\delta =1/{\rm{\Gamma }}(1-\beta \cos {\theta }_{{\rm{v}}})$ (see, e.g., Ghisellini et al. 2006). To this end, we sample the Γ–E_iso correlation reported by Ghirlanda et al. (2018), to assign Γ values to the simulated events.

The GRB prompt emission spectrum is described with the empirical Band function (Band et al. 1993), consisting of two smoothly joined power laws. To each simulated GRB we assign a low- and high-energy spectral index, extracted from the Gaussian distributions ${\mathfrak{G}}[\alpha | \mu ,\sigma ]$ and ${\mathfrak{G}}[\beta | \mu ,\sigma ]$, respectively, whose parameters (e.g., Goldstein et al. 2012) are given in Table 1.

The population model has 12 free parameters (Table 1): [ρ₀, p_z,1, p_z,2, p_z,3] describe the GRB cosmic formation rate (Equation (1)) and [a₁, a₂, L_iso,0, δ] define the luminosity function (Equation (2)); the remaining four parameters (μ, σ) and (K, S) describe the stretching factor of the duration distribution and the E_iso–θ_jet correlation, respectively. All the other parameters reported in Table 1 are fixed to their presently known values and are derived from the corresponding references given above.

The 12 free parameters are constrained by matching a set of observed distributions C_J (Section 3) and by maximizing a binned likelihood. The probability that the model distribution has m_i objects in the i bin of the distribution, given that there are r_i real events in this bin, follows a Poisson distribution. The total log-likelihood can be expressed as

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\sum _{{C}_{j}}\sum _{i=1}^{n}{r}_{i}\mathrm{ln}({m}_{i})-{m}_{i}-\mathrm{ln}({r}_{i}!),\end{eqnarray} \tag{ 3 }$$

where the internal sum is performed over n bins of each distribution of the constraint C_J (with J = 1,...,14). The observational constraint distributions C_J (described in Section 3) are built from the observed properties of the LGRB samples detected by different instruments, and are represented by the different symbols in the panels of Figure 1.

Zoom In Zoom Out Reset image size

The posterior distributions of the 12 free parameters are derived through a Markov Chain Monte Carlo (MCMC) method, which samples from wide uniform prior distributions of the free parameters: p_z,1 ∈ [1, 6], p_z,2 ∈ [0.1, 6], p_z,3 ∈ [1, 9], δ ∈ [0, 5], a₁ ∈ [0.5, 3], a₂ ∈ [2, 5], logL_iso,b ∈ [50, 53], μ_f ∈ [1.5, 10], σ_f ∈ [0.01, 6], K_E−θ ∈ [9, 10.5], and S_E−θ ∈ [ −0.25, −0.1]. We implement the Goodman & Weare (2010) MCMC algorithm. We work with 200 walkers running over 200 steps, and we adopt the parallel stretch move method (Foreman-Mackey et al. 2013) to reduce the computation time.

The free parameters of the GRB population model are constrained by reproducing a set of distributions of the observables provided by samples of GRBs detected by different missions. The numerous samples of LGRBs detected by the Gamma Burst Monitor (GBM) on board the Fermi satellite and the Burst And Transient Source Experiment (BATSE) on board the Compton Gamma Ray Observatory (CGRO) provide the peak flux (F_P ) and fluence (F) distributions (panels (a) and (c) of Figure 1), as well as the observer-frame peak energy (E_p/1 + z) and duration (T₉₀) (panels (b) and (d) of Figure 1). In order to minimize sample incompleteness, which is particularly severe at the faint end of the peak flux distribution, ⁴ we consider LGRBs with F_P ≥ 4 ph cm⁻² s⁻¹, integrated over the 10 keV–1 MeV energy range for Fermi/GBM and the 20 keV–2 MeV range for GCRO/BATSE. This selection results in 904 GRBs (out of 1906) for Fermi/GBM (the red symbols in Figure 1) and 370 (out of 1529) for CGRO/BATSE (the blue symbols in Figure 1). The inset in Figure 1(a) shows the peak flux (50–300 keV) of 2810 CGRO/BATSE bursts in the Stern et al. (2001) catalog, with F_p ≥ 0.3 ph cm⁻² s⁻¹.

For the redshift distribution (panel (e)), we consider two well-selected Swift/BAT flux-limited samples, in which large majorities of the bursts have a redshift measurement. The BAT6 sample (Salvaterra et al. 2012; Pescalli et al. 2016) consists of LGRBs with a 15–150 keV peak flux ≥2.6 ph cm⁻² s⁻¹. The SHOALS sample (Perley et al. 2016) comprises GRBs with a 15–150 keV integrated fluence ≳10⁻⁶ erg cm⁻². Both these samples have 80%–90% redshift measurements, either from afterglow or host spectroscopy. The BAT6 sample also provides the constraints on the E_iso and L_iso distributions (panel (f)). Finally, we consider the jet opening angle distribution (panel (g)) of the GRBs collected in Ghirlanda et al. (2007).

Further, we also verified the consistency of our best model (Figure 5 in Appendix A.1) with the distributions of the peak flux, fluence, and duration (T₉₀) of the samples of LGRBs observed by Swift/BAT ⁵ and the Gamma Ray Burst Monitor (GRBM) on board the Beppo/SAX satellite (Frontera et al. 2008).

The model that best matches the set of observational constraints is shown by the solid histograms in Figure 1 (and in Figure 5). The shaded/dashed regions around the model line represent the 68% model uncertainty. This is computed by resampling the joint posterior distributions of the free parameters (shown in Figure 2) with 10⁵ samples, then computing the corresponding model for each. For each fixed value of an independent variable, the shaded region shows the 16th to 84th percentile of the resulting model. The posterior distributions of the best model are shown by the corner plot in Figure 2, and their central values and 1σ uncertainty are reported in Table 1. All the parameter values are well constrained with, on average, ∼10% accuracy.

Zoom In Zoom Out Reset image size

The luminosity function has a faint-end slope ∼ −1 and a bright-end slope ∼ −2.2, with a break in luminosity at z = 0 L_b,0 = 10⁵² erg s⁻¹. These values are consistent with those reported by previous studies (Wanderman & Piran 2010; Salvaterra et al. 2012). We find a mild evolution of this break, which increases with redshift as L_b = L_b,0(1 + z)^0.64.

The cosmic rate of the LGRBs, ρ(z), is shown in Figure 3(a) (the solid red line). The shaded orange region represents the 1σ uncertainty, as obtained by resampling the posterior distributions of the free parameters. ρ(z) increases with redshift ∝(1 + z)^3.2, it peaks at z ∼ 3, and it declines at larger redshifts ∝(1 + z)⁻³. These slopes are steeper and similar, respectively, when compared to those of the CSFR ⁶ of Madau & Dickinson (2014) (the solid blue line in Figure 3(a)).

Zoom In Zoom Out Reset image size

The local LGRB formation rate is ${\rho }_{0}\,=\,{79}_{-33}^{+57}$ Gpc⁻³ yr⁻¹ (68% confidence interval). This is derived by normalizing the model population: we consider the ∼170 ± 13 bright GRBs with peak fluxes in the 15–150 keV energy range ≥2.6 ph cm² s⁻¹ detected by Swift/BAT over 8 yr. We assume for Swift/BAT an average duty cycle of 0.75 and a fully coded field of view of 1.4 sr. We have verified that our results are unchanged if we renormalize the population to the bright GRB sample detected by Fermi/GBM, namely 902 LGRBs with P(10–1000) ≥ 5 ph cm² s⁻¹, detected over 11.8 yr within its 8.8 sr field of view, also accounting for a 0.55 duty cycle. Since our population accounts for the random orientation of GRB jets in the sky, the derived local rate represents the intrinsic one.

We also performed the same analysis under the assumption of a Schechter luminosity function, and the results are reported in Appendix A.2. We find that the shape of the LGRB rate at low redshifts is unchanged (Figure 6). The Schechter scenario produces a slightly shallower decrease in the GRB formation rate at high redshifts. However, this is consistent with the CSFR, given the uncertainties on high-redshift measurements.

The LGRB efficiency as a function of redshift, defined as the ratio between the GRB rate and the CSFR, is shown in Figure 3(b) by the solid red line. Its 1σ uncertainty (the shaded violet region in Figure 3) is computed by resampling the posterior distributions of the free parameters shown in Figure 2 10⁵ times. For the CSFR, we have also considered the function of Madau & Dickinson (2014; the solid blue line in Figure 3(a)). In the local universe, this ratio is ∼5 ×10⁻⁵ M_⊙ ⁻¹, and it increases by up to a factor of 10 at z ∼ 3. In Figure 3(b), we show with the blue and red symbols, respectively, the fraction of star formation leading to Type Ibc and Type Ic broad-line (BL) supernovae (Perley et al. 2020), assuming the core-collapse supernova rate derived by Strolger et al. (2015). We find that, in the local universe, ${1.3}_{-0.6}^{+1.0} \%$ of Type Ic-BL supernovae produce a jet and the associated GRB event. This fraction increases with redshift up to ∼7% at z > 3.

The general trend of the GRB cosmic rate at low redshifts shows a decrease in the efficiency of GRBs being produced with decreasing redshift (Figure 3). This result argues against the claims for an increase in the LGRB efficiency toward lower redshifts (the so-called low-redshift excess—Yu et al. 2015; Petrosian et al. 2015; Tsvetkova et al. 2017; Lloyd-Ronning et al. 2019; but see also Pescalli et al. 2016; Le et al. 2020; Bryant et al. 2021) obtained through nonparametric methods. The low-redshift GRB uptick could be partially absorbed if the jet opening angle evolves with redshift (Lloyd-Ronning et al.2020).

The GRB efficiency with redshift evolution can be interpreted as being due to the increase in metallicity with cosmic time, which prevents a larger fraction of massive stars from producing GRBs (Yoon et al. 2006). In order to verify this possibility, we compute the fraction of cosmic star formation in galaxies with average metallicity below a given threshold, following the method of Robertson & Ellis (2012). We combine (Appendix A.3) the mass–metallicity relation and its redshift evolution (Maiolino et al. 2008) with the galaxy mass function (McLeod et al. 2021) and the galaxy main sequence (Tomczak et al. 2014). We find that a metallicity threshold $12+\mathrm{log}({\rm{O}}/{\rm{H}})\leqslant$ 8.6, represented by the the dashed blue lines in Figure 3, accounts for the derived GRB efficiency evolution. The metallicity thresholds we derive are consistent with the observed GRB host galaxy metallicity (Vergani et al. 2015; Palmerio et al. 2019). To further verify this consistency, we have computed the stellar masses of galaxies with metallicities below the threshold value found, and compared them with the observed masses of a complete sample of GRB hosts (Vergani et al. 2015; Palmerio et al. 2019; the black symbols in the left panel of Figure 4). The maximum and average stellar masses (the green and blue curves in the left panel of Figure 4) are shown for $12+\mathrm{log}({\rm{O}}/{\rm{H}})\,\leqslant$ 8.6 by the dashed lines and they are consistent with the measured masses of the hosts and their redshift evolutions.

Zoom In Zoom Out Reset image size

Previous LGRB population studies (e.g., Firmani et al. 2004; Salvaterra et al. 2012; Palmerio & Daigne 2021) have showed that a pure density evolution (PDE) of the GRB formation rate, where the luminosity function is invariant with redshift, and a pure luminosity evolution (PLE), where the GRB formation rate is proportional to the CSFR, cannot be distinguished on the basis of the GRB peak flux and luminosity–redshift distributions. Our best model is intermediate between these two extreme scenarios, and, indeed, we find both the GRB formation rate evolving with redshift (Figure 3(b)) and a mild evolution of the characteristic GRB luminosity ∝(1 + z)^0.64 (Table 1). We exclude at more than 5σ significance the extreme cases of PLE and PDE. Moreover, the PLE is also excluded by the host galaxy mass distribution (see the red line in Figure 4). The finding of the evolution of the LGRB characteristic luminosity ∝(1 + z)^0.64 could hide the changes in some properties of the progenitors with redshift, possibly ruled by their metallicity (Yoon et al. 2006), which could induce a flattening of the initial mass function (Fryer et al. 2022).

This possibility can be tested by the increasing fraction of GRBs detected at high redshift, as envisaged by the next-generation missions optimized for this purpose (e.g., the Transient and High Energy Sky and Early Universe Surveyor, or THESEUS—Amati et al. 2021; or the Gamow Explorer—White et al. 2020).

Our population model, which implements the jet opening angles of GRBs and accounts for the observed properties as a function of the viewing angle, allows us to estimate the fraction of GRBs detected by Fermi/GBM and Swift/BAT that are observed off-axis, i.e., θ_v > θ_jet. The right panel of Figure 4 shows the expected cumulative rate as a function of redshift for Fermi/GBM (blue filled stripes) and Swift/BAT (orange filled stripes). The hatched curves show the detection rates of off-axis events according to our population. We estimate that ∼10% of the LGRBs detected by Fermi and Swift at z < 0.5, corresponding to a rate of 0.2–0.7 yr⁻¹, should be observed off-axis. This fraction reduces to less than 3% at z < 2. In this work, we have adopted a uniform jet structure consistent with that suggested by study of the luminosity function of LGRBs (Pescalli et al. 2015). Therefore, our results should be considered as a conservative lower limit, since a universal structured jet should produce a larger fraction (e.g., Salafia et al. 2015) of detected off-axis events.

We acknowledge support from ASI-INAF agreement No. 2018-29-HH-0. We acknowledge PRIN-MIUR 2017 (grant 20179ZF5KS) and the "Figaro" premiale project (1.05.06.13). We thank S. Vergani, G. Ghisellini, and O. S. Salafia for constructive discussions. We acknowledge L. Amati and N. White, and the respective THESEUS and Gamow Explorer collaborations, for stimulating interactions over the past few years.

A.1. Consistency Checks

In order to verify the consistency of our model with the GRB samples detected by other instruments, we considered the full Swift/BAT and BeppoSAX/GRBM (Frontera et al. 2008) samples. We applied a peak flux cut, extracting 536 Swift/BAT GRBs with a 15–150 keV peak flux F_p ≥ 0.5 phot cm⁻² s⁻¹. Their observer-frame distributions are shown by the green asterisks in the top panels of Figure 5. For Beppo/SAX, we considered the 118 GRBs with a peak flux, integrated in the 40–700 keV energy range, F_p ≥ 10⁻⁶ erg cm⁻² s⁻¹. The Beppo/SAX distributions are shown by the pink asterisks in the bottom panels of Figure 5. The population model is shown by the solid histogram and shaded regions.

Zoom In Zoom Out Reset image size

A.2. Schechter Luminosity Function

We considered the possibility that the GRB luminosity function is described by a Schechter function (Schechter 1976), namely a power law with a high-luminosity exponential cutoff. This model has two free parameters, the slope of the power law α and the cutoff luminosity L_c . We fix the minimum luminosity of this function to 10⁴⁷ erg s⁻¹. With this model, we can reproduce all the observational constraints described in Section 3. We obtain α = 1.15 ± 0.15 and $\mathrm{log}({L}_{c})=53.21\pm 0.28$. The resulting GRB formation rate is shown in Figure 6 by the black line (and its 1σ uncertainty by the shaded region). For comparison, the CSFRs of Madau & Fragos (2017; blue line) and Madau & Dickinson (2014; green line) are shown. The modified CSFRs obtained with metallicity threshold cuts at 12+log(O/H) ≤ 8.6 and 8.4 are shown by the dashed and dotted–dashed lines, respectively. The GRB formation rate obtained under the Schechter assumption is consistent within its uncertainties with that obtained with a broken-power-law luminosity function, although it is slightly shallower at high redshifts.

Zoom In Zoom Out Reset image size

A.3. Metallicity Bias Model

Following the formalism of Robertson & Ellis (2012), we estimate the fraction of star formation occurring in galaxies with metallicity below a given threshold as

$$\begin{eqnarray}&&f(z)=\displaystyle \frac{{\int }_{{M}_{\star ,\min }}^{{M}_{\star ,\mathrm{th}}(z)}{R}_{\star }(M,z)\cdot {\rm{\Upsilon }}(M,z){dM}}{{\int }_{{M}_{\star ,\min }}^{{M}_{\star ,\max }}{R}_{\star }(M,z)\cdot {\rm{\Upsilon }}(M,z){dM}},\end{eqnarray} \tag{ A1 }$$

where R_⋆(M, z) is the star formation rate–stellar mass relation (M is the galaxy stellar mass in units of M_⊙) and ϒ(M, z) is the galaxy mass function. Both these functions evolve with the redshift z. For R_⋆(M, z), we assume the relation of Tomczak et al. (2014), and for ϒ(M, z) we consider the recently reported double Schechter galaxy mass function of McLeod et al. (2021). M_⋆,th(z) is the maximum mass corresponding the the metallicity threshold,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{M}_{\star ,\mathrm{th}}(z) & = & {\left[\displaystyle \frac{K(z)-\left[12+\mathrm{log}({\rm{O}}/{\rm{H}})\right]}{0.0864}\right]}^{1/2}\\ & & +\mathrm{log}{M}_{0}(z)-0.04,\end{array}\end{eqnarray} \tag{ A2 }$$

obtained by inverting Equation (2) of Maiolino et al. (2008). The values of K(z) and M₀(z) are interpolated from those reported in Table 5 of their paper. We also corrected (through the 0.04 term) for the different assumption of a Chabrier initial mass function in the host galaxy mass estimates shown in Figure 4.