The impact of the FMR and starburst galaxies on the (low metallicity) cosmic star formation history Open Access

ABSTRACT

1 INTRODUCTION

The observed populations of stars, their remnants and related transients consist of objects formed at different times and with different metallicities. To model those populations and to correctly interpret observations, it is necessary to know the metallicity dependent star formation history (SFH) of the observed galaxy, galaxies in the probed volume, or even of the entire Universe. Knowledge of the latter – which is the subject of interest of this study – becomes increasingly important in the era of gravitational wave (GW) astrophysics. That is because the time between the formation of the progenitor stars and merger of stellar black holes or neutron stars observed in GW can be comparable to the age of the Universe (e.g. Belczynski et al. 2016). Moreover, the efficiency of formation of merging binaries may show a strong metallicity dependence. In particular, it has been suggested that double black hole mergers may form much more efficiently in low metallicity environments (≲0.1 solar metallicity; e.g. Giacobbo, Mapelli & Spera 2018; Klencki et al. 2018). This makes the modelled properties of the population of such systems particularly sensitive to the assumed distribution of the cosmic star formation rate density (SFRD) at different metallicities and redshifts, f_SFR(Z,z) (e.g. Chruslinska, Nelemans & Belczynski 2019; Neijssel et al. 2019; Broekgaarden et al. 2021) and requires knowledge of f_SFR(Z,z) even beyond the peak of the cosmic SFH. To confront the model predictions with observations and draw correct conclusions from such a comparison, it is necessary to take into account the f_SFR(Z,z) uncertainty, which may be substantial – especially at high redshifts (Chruslinska & Nelemans 2019; Chruślińska et al. 2020; Boco et al. 2021). Given the many uncertain pieces of information that need to be combined to estimate f_SFR(Z,z), an (observation-based) determination of this distribution presents a challenge in itself. In this study, we take a closer look at two of those pieces: the empirical correlation between the star formation rate and metallicity of star-forming galaxies (the so-called fundamental metallicity relation, e.g. Ellison et al. 2008; Mannucci et al. 2010) and the contribution of starburst galaxies. We build on the observation-based f_SFR(Z,z) model from Chruslinska & Nelemans (2019) (briefly introduced in Section 2) and expand the discussion presented in Boco et al. (2021), aiming to evaluate the uncertainty (or – find realistic, observationally allowed extremes) of the metallicity dependent cosmic SFH in view of those factors. Where appropriate we adopt a standard flat cosmology with Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹, assume a (universal) Kroupa (2001) initial mass function (IMF) and solar metallicity of 12 + log₁₀(O/H)_⊙ = Z_O/H⊙ = 8.83 and Z_⊙ = 0.017 (Grevesse & Sauval 1998).

2 AN OBSERVATION BASED f_SFR(Z, z) DETERMINATION

We construct f_SFR(Z,z) building on the framework detailed in Chruslinska & Nelemans (2019). In essence, their method is based on a combination of three key ingredients: the distribution of galaxy stellar masses (galaxy stellar mass function – GSMF), and the distributions describing the star formation rates (SFR) and metallicities of galaxies at fixed stellar mass M_*. All distributions are redshift dependent. The GSMF is used to obtain the number density of galaxies of different masses. The contribution of galaxies of different masses to the total SFRD can then be obtained by weighing the number density of galaxies in a given mass range by their SFR. Finally, f_SFR(Z,z) is obtained by assigning a metallicity to each M_*.

The SFR and metallicity distributions are modelled as lognormal distributions centred on the empirical star formation–mass relation (SFMR) and the mass – (gas-phase) metallicity relation (MZR), respectively, with dispersions σ_SFMR and σ_MZR (representing the intrinsic scatter around the two relations). Note that the observational gas-phase metallicity estimates probe the oxygen abundance (Z_O/H = 12 + log₁₀(O/H)¹), and that is the metallicity measure used in our study. Additional scatter in Z_O/H is introduced to model the spread in the metallicity at which the stars are forming within the galaxies. To determine all the necessary ingredients, the authors assemble a compilation of observational results describing the MZR, SFMR and GSMF, combined over a wide range of redshifts (z) and M_*. Several variations of the base relations (SFMR, MZR, GSMF) are explored in order to discuss the impact of the uncertain absolute metallicity scale (steming from the differences between the estimates obtained with different metallicity determination methods), the shape of the high-mass end of the SFMR and the redshift evolution of the low-mass end of the GSMF (see Section 2 in Chruslinska & Nelemans 2019 and references therein).

In this study, we explore an alternative method to obtain the redshift-dependent metallicity distribution of star-forming galaxies, based on the fundamental metallicity relation. We then modify the SFR distribution to better account for the contribution of star-forming galaxies that are strong outliers to the general SFMR (starbursts). Our motivation for those modifications and the details of our approach are laid out in Sections 3 and 4, respectively.

3 THE FUNDAMENTAL METALLICITY RELATION

The fundamental metallicity relation (FMR) is a three parameter dependence linking M_*, SFR, and gas-phase metallicity of star-forming galaxies (Ellison et al. 2008; Mannucci et al. 2010). The relation implies the mass–metallicity correlation, as observed in the MZR, but it also introduces an anticorrelation between metallicty and SFR, such that galaxies of the same stellar mass showing higher than average SFR also have lower metallicities (i.e. the galaxy’s offset from the average mass–metallicity relation – MZR – is anticorrelated with its offset from the average SFR–mass relation).

The existence of such SFR–metallicity anticorrelation has been reported in numerous observational studies (e.g. Lara-López et al. 2010; Mannucci, Salvaterra & Campisi 2011; Yates, Kauffmann & Guo 2012; Andrews & Martini 2013; Bothwell et al. 2013; Salim et al. 2014; Zahid et al. 2014b; Yabe et al. 2015; Hunt et al. 2016; Sanders et al. 2018; Cresci, Mannucci & Curti 2019; Curti et al. 2020; Sanders et al. 2021) and established up to redshift z ∼ 3.5.

Such correlation is also expected from theoretical models of galaxy evolution and is thought to reflect the changes in the SFR-fuelling gas fraction present within the galaxy (e.g. Ellison et al. 2008; Davé, Finlator & Oppenheimer 2011; Yates et al. 2012; De Lucia et al. 2020). Lower SFR for a given M_* then implies that the galaxy has already used up most of its cold gas reservoir (fuelling its SFR in the past), and therefore (in the absence of strong outflows) its interstellar medium is more metal rich than that of the galaxy of the same mass that is still highly star forming. At low SFR there are also fewer supernovae that can remove (especially at low M_*) metal-rich gas from the galaxy. In turn, inflowing metal-poor material lowers the metallicity and provides additional fuel for the SFR, increasing the latter. As long as the time-scales on which the SFR and metallicity evolve are of the same order, the anticorrelation between the galaxy offsets from the MZR and SFMR is expected to hold (e.g. Torrey et al. 2018). Those time-scales are similar for the feedback/SFH implementations used in the large-scale cosmological simulations such as EAGLE and Illustris-TNG, which warrants the existence of the FMR up to high redshifts in the simulations (Lagos et al. 2016; De Rossi et al. 2017; Torrey et al. 2018). However, Torrey et al. (2018) point out that the strength of the correlation – especially for low-mass galaxies at high redshift – may be reduced if models with particularly strong and/or bursty feedback are used. Lagos et al. (2016) find that the shape of the FMR in the EAGLE simulations depends on the adopted model of star formation. Therefore, observational confirmation of the existence or breakdown of the FMR at high redshifts will help to discriminate between the different feedback and star formation prescriptions.

Observationally, the FMR is found to show little to no evolution with redshift up to z ∼ 3 within the uncertainty of the current data and range of galaxy properties probed (see Cresci et al. 2019 and Sanders et al. 2021 for recent discussion). Since both the mass and SFR distributions of star-forming galaxies change over time, the apparent lack of FMR evolution implies that galaxies at different redshifts probe different parts of the locally established relation. At fixed M_*, a decrease in the average galaxy metallicity as a function of redshift is still expected, as the typical SFR is higher at earlier cosmic times.

At z ≳ 2, the rate of decrease in metallicity implied by a non-evolving FMR is much weaker than that reported by various studies discussing the redshift evolution of the MZR (e.g. Maiolino et al. 2008; Mannucci et al. 2009). Similarly to the FMR, the MZR (and its evolution with redshift) is virtually unconstrained at redshifts z ≳ 3.5 (Maiolino & Mannucci 2019). This uncertainty in the rate of metallicity evolution is one of the key factors affecting the f_SFR(Z,z). In particular, different high redshift extrapolations of empirical relations used in the literature lead to drastically different conclusions about the birth metallicities of stars forming beyond the peak of the cosmic SFH (e.g. Chruslinska et al. 2019; Chruslinska & Nelemans 2019; Boco et al. 2021).

Recent observational studies find support for rapid early metal enrichment in high-redshift galaxies, pointing towards a weak MZR evolution that is more in line with the apparently invariant FMR (see section 3 in Boco et al. 2021 and references therein for a recent discussion). Such weaker metallicity evolution is also in agreement with the results of cosmological simulations and semi analytical models of galaxy evolution (e.g. Yates et al. 2012; Torrey et al. 2019). Furthermore, Sanders et al. (2021) show that when the differences in the properties of the ionized gas within H ii regions (responsible for the emission lines used to estimate the metallicity; in particular lower iron to oxygen ratio for the same O/H found in high-redshift galaxies) of galaxies at different redshifts are taken into account, the inferred MZR evolution is milder (i.e. metallicities of high-redshift galaxies are biased low if those differences are neglected). In fact, the MZR evolution found by Sanders et al. (2021) is consistent with the non- or weakly evolving FMR within the redshift and mass range probed in their study. The existence of an FMR also means that if the galaxy sample is biased towards high SFRs, the inferred average metallicity is underestimated. Such biases can be expected in high redshift and low stellar mass galaxy samples, affecting the MZR shape (the low mass end slope) and the rate of evolution with redshift (a decrease in the normalization).

In light of the above discussion, we conclude that the assumptions about the metallicity evolution made in Chruslinska & Nelemans (2019), which rely on the MZR obtained by Mannucci et al. (2009) and extrapolate its redshift evolution based on their two highest redshift bins, are likely to overestimate the rate of metallicity decrease at z ≳ 2 and require revision: either assuming a weaker MZR evolution, or using the FMR to assign metallicity to galaxies. The use of the FMR instead of an MZR in f_SFR(Z,z) determination appears advantageous for two reasons: (i) the redshift invariance of the FMR allows to circumvent the problem of the uncertain MZR evolution with redshift - the metallicity distribution at each z is then defined by the local FMR and the SFMR and GSMF (whose z dependence is generally better constrained than that of the MZR). We stress that with the current data there is no guarantee that at z ≳ 3 the FRM remains (close to) redshift invariant. By assuming a non-evolving FMR throughout the entire cosmic history we explore another extreme assumption (with respect to Chruslinska & Nelemans 2019) about the rate of metallicity evolution at high redshift. (ii) if all star-forming galaxies follow the FMR (at present, there is no clear evidence to the contrary), it can be used to describe metallicity of galaxies belonging to populations for which there are little examples of metallicity determinations (e.g. strong outliers of the star-forming main-sequence that are not described by the MZR, such as starburst galaxies, see Section 4).

However, the major problem with this approach is that there is no agreement on the exact form of this three parameter mass–SFR-metallicity dependence. The f_SFR(Z,z) derived with the FMR will necessarily strongly depend on the choice of this relation, just as it strongly depends on the choice of MZR and its extrapolated evolution. While the discussion of the systematic effects on the f_SFR(Z,z) introduced by the choice of a particular form of the MZR is relatively straightforward (Chruslinska & Nelemans 2019) and to a large extent boils down to the discussion of differences caused by the use of different metallicity determination methods (which are well documented in the literature, e.g. Kewley & Ellison 2008, Maiolino & Mannucci 2019), an analogous discussion for the FMR is not. The observationally inferred FMR is known to non-trivially depend on the metallicity determination method, the SFR determination method and the galaxy sample selection criteria (e.g. Yates et al. 2012; Hunt et al. 2016; Kashino et al. 2016; Telford et al. 2016; Cresci et al. 2019). We note that it is not clear how the FMR depends on each of those factors. Therefore, one needs to be cautious when using different literature results concerning the properties of galaxy populations (e.g. the SFMR) in combination with a particular FMR estimate – if those results are based on different techniques, this would represent an internally inconsistent approach.

Instead of using a particular example from the literature, in our analysis we resort to a phenomenological description of the FMR. We aim to capture the robust observational/theoretical features of the mass–SFR-metallicity dependence and describe the variations in the shape of this dependence caused by different choices of the local MZR and SFMR (that reflect the differences in the metallicity and SFR determination methods used in the literature). This allows us to ensure the consistency of the method and discuss the realistic extremes of the f_SFR(Z,z).

3.1 Constructing the FMR from z ∼ 0 scaling relations

To construct our model FMR, we describe its 2D projection on the Z_O/H–log₁₀(M_*) plane. This projection is the most commonly shown in observational studies. We refer to it as Z_O/H(SFR, M_*) (see Fig. 1 for the illustration).

Figure 1.

Sketch illustrating the local MZR (top left), SFMR (bottom left), and the three parameter mass–metallicity-SFR dependence Z_O/H(SFR, M_*) (right), indicating the key parameters of those relations and the observed anticorrelation between the galaxy’s offsets from the MZR (δ_MZR) and SFMR (δ_SFMR). We consider different assumptions about the MZR and SFMR relations (i.e. normalization |$\rm \mathit{ Z}_{O/H; MZR0}$| and slope a_MZR of the MZR, high-mass end of the SFMR a_SFMR2) to cover the range of possibilities present in the literature, as indicated with the thick grey lines in the left-hand panels. In Section 3.1, we relate the parameters of the Z_O/H(SFR, M_*) to the parameters of the local 2D relations.

Open in new tabDownload slide

3.1.1 The asymptotic metallicity

A clear feature present among the observational Z_O/H(SFR, M_*) determinations is the high mass flattening. We note that some authors report a reversal of the (anticorrelation at high M_* rather than a simple flattening, e.g. Yates et al. 2012; Kashino et al. 2016). Telford et al. (2016) point out that certain sample selection criteria (e.g. applying signal-to-noise ratio cuts on oxygen lines used to estimate the metallicity) may lead to biases against the massive, low SFR galaxies and therefore affect the inferred low SFR and high mass part of the FMR. They argue for such biases as the cause of the reversal of the correlation at high masses as seen in Kashino et al. (2016). The effects of dust (and the applied dust SFR corrections) may also induce biases against massive, metal-rich galaxies (Telford et al. 2016).

Change in the relation at high masses is also supported theoretically, and could be attributed to the increased importance of AGN feedback that can rapidly influence the SFR while the metallicity continues to evolve much more gradually (e.g. De Rossi et al. 2017; Torrey et al. 2018). We therefore assume that the high mass flattening is a robust feature of the Z_O/H(SFR, M_*) and include it within our model. Observationally, the value of the asymptotic metallicity Z_{O/H; 0} at which the relation flattens appears to be roughly independent of the SFR and coincide with that of asymptotic metallicity of the average z ∼ 0 MZR (Z_{O/H; MZR0}). As such, it is affected by the choice of the metallicity determination method, which leads to systematic offsets in the normalization of the MZR (compare examples shown in Telford et al. 2016 and in fig. 2 in Cresci et al. 2019). We fix Z_{O/H; 0} = Z_{O/H; MZR0}, where Z_{O/H; MZR0} is defined by the choice of the z ∼ 0 MZR.

Figure 2.

Dependence of the parameters γ [Z_O/H(SFR, M_*) slope at M_* << M_{0; SFR}, see equation (1)] and α [slope of the log₁₀(M_{0; SFR}) dependence on log₁₀(SFR), see equation (6)] on the slopes of the z ∼ 0 MZR and SFMR and the parameter ∇_FMR0. Hatched areas in the top panels span between the values obtained for a_MZR ∼ 0.3 (orange lines/areas) and 0.6 (blue lines/ares) and indicate the variation caused by the different MZR slopes considered in this study. Grey areas in the top panels and in the bottom right panel span between a_SFMR = 0.7−1 to illustrate the variation caused by the change in the SFMR slope. The hatched green area in the bottom left panel shows the variation caused by varying ∇_FMR0 in the range 0.17–0.3, assuming a_SFMR = 0.83 (our default low/intermediate mass slope of the SFMR). The dashed green line falling within that range indicates our fiducial choice of ∇_FMR0 = 0.27. Orange (blue) ranges in the bottom right panel show the range of α obtained with a_MZR ∼ 0.3 (a_MZR ∼ 0.6) and ∇_FMR0 in the range 0.17–0.3 for each a_SFMR. The orange square shows α and a_MZR as found by Curti et al. (2020). The brown diamond shows α/γ calculated with a_MZR, a_SFMR, and ∇_FMR0 from Sanders et al. (2021) (note that it overlaps with the estimate from Curti et al. 2020 in the bottom left panel).

Open in new tabDownload slide

3.1.2 The slope at low masses and high SFR

3.1.3 Turnover mass as a function of SFR

Observational studies suggest that |$M_{0;\rm SFR}$| increases with SFR [i.e. for higher SFR, Z_O/H(SFR, M_*) flattens at higher masses]. The steepness of this dependence governs the spacing between the Z_O/H(M_*,SFR) curves obtained for different SFR (see Fig. 1). Examples shown in Telford et al. (2016) and in fig. 2 of Cresci et al. (2019) suggest that both the metallicity and SFR determination methods have impact on the distance between the different SFR lines (and so on M_0;SFR). Curti et al. (2020) find that a linear dependence (log₁₀(M_0;SFR) ∝ log₁₀(SFR)) can well describe the trend seen for their z ∼ 0 galaxy sample.

3.1.4 The value of β

β describes the bending of the Z_O/H(SFR, M_*) relation. The higher the value of β, the sharper the transition from the linear to flat regime and the smaller the range of log₁₀(M_*) over which the transition happens (see Fig. A1 for illustration). β could in principle depend on the SFR. Curti et al. (2020) let β as a free parameter when fitting their Z_O/H(SFR, M_*) relation for different SFR bins. While their best-fitting relations give different β values in different SFR bins, the authors note that there is no clear dependence on SFR. Ultimately, they quote a single value of β ∼ 2 to describe the FMR.

In principle, we can numerically solve for β using the condition that galaxies with zero offset from the z ∼ 0 SFMR lie on the z ∼ 0 MZR. This condition allows us to properly recover the z ∼ 0 MZR when sampling the galaxy properties from the z ∼ 0 GSMF and SFMR and assigning their metallicities with Z_O/H(SFR, M_*). This condition is only applicable in the SFR range that is covered by the z ∼ 0 SFMR (≳10⁻³ M_⊙ yr⁻¹ and ≲10² M_⊙ yr⁻¹). Numerous solutions are allowed – especially at low SFR ≲ 10⁻¹ M_⊙ yr⁻¹, where the intersection with the z ∼ 0 MZR is reached at log(M_*) <log(⁠|$M_{0;\rm SFR}$|⁠) i.e. where Z_O/H(SFR, M_*) is virtually insensitive to the choice of β. We use the numerical solutions as a guide and adopt a simplified β dependence on log₁₀(SFR), verifying that Z_O/H(SFR, M_*) for galaxies at z ∼ 0 SFMR leads to Z_O/H that agree with the corresponding values from the z ∼ 0 MZR to within 0.01 dex (note that this is smaller than the typically found residual scatter about the Z_O/H(SFR, M_*) of ∼0.05 dex). Where necessary, we extrapolate the obtained dependence to lower and higher SFR values than probed by z ∼ 0 SFMR. The resulting dependence is shown in Fig. A2 in Appendix  A.

3.1.5 The strength of the SFR–metallicity correlation

Observationally, ∇_FMR appears to be a function of M_* and SFR. Such dependence is evident in the analysis presented by Salim et al. (2014) (see Table 1 therein, where κ corresponds to our ∇_FMR), where the authors compare the strength of the correlation for galaxies split in different mass and SFR bins. This is one of the key observed characteristics of the FMR that we aim to reproduce within our description.

The fact that Z_{O/H; 0} appears to be roughly independent of the SFR means that the observed Z_O/H(SFR, M_*) flattening does not simply reflect the presence of the high mass flattening of the MZR (and potentially SFMR), but also indicates a weakening/disappearance of the SFR–Z_O/H correlation in the high M_* regime (i.e. small/zero ∇_FMR). Note that, by construction (fixing the value of Z_{O/H; 0}), we recover such high M_* weakening of the correlation in our method. A similar behaviour is also seen in the Z_O/H – log₁₀(SFR) projection of the FMR, where at fixed log₁₀(M_*) the flattening appears at relatively low SFR values (compared to the SFMR). In other words, the correlation between the SFR and metallicity is found to weaken/disappear at low SFR/low specific SFRs (sSFR = SFR/M_* ≲ 10⁻¹⁰–10^−10.5 yr⁻¹). When projected on to the Z_O/H–log₁₀(M_*) plane, this means that the correlation is weaker ‘above’ than ’below’ the average z ∼ 0 MZR. That is also the case within our description: for the same absolute value of δ_SFMR, ∇_FMR is smaller above than below the z ∼ 0 MZR due to the SFR-dependent location of the flattening and – in certain model variations – the β dependence on SFR.

While observational studies seem to agree on the existence of those high M_* and low SFR/low sSFR weaker correlation regimes, there is no agreement on the precise region of the FMR [in terms of the Z_O/H, log₁₀(M_*), and log₁₀(SFR) values] in which they appear. Within our description we can qualitatively reproduce those trends. The transition between the strong correlation regime and the part of the FMR where the correlation weakens happens at different log₁₀(M_*) and sSFR depending on the parameters of the z ∼ 0 MZR and SFMR. We stress that we only assume that ∇_FMR is fixed (i.e. ∇_FMR = ∇_FMR0, independent of M_* and SFR) in the low mass and high SFR part of the FMR – as guided by the observed characteristics of Z_O/H(SFR, M_*) (but see Appendix  C, where we further relax this assumption and introduce an SFR-dependent ∇_FMR in this region of the FMR). We do not explicitly assume any value or dependence for ∇_FMR in the remaining part of the relation.

∇_FMR0 is the only parameter used in our description of the Z_O/H(SFR, M_*) that is not defined by the choice of the local MZR and SFMR, as a potential dependence of this parameter on the metalicity/SFR derivation method is unclear and only few determinations of ∇_FMR are given in the literature. Using the H α based SFR and metallicity derivation method of Mannucci et al. (2010) and Salim et al. (2014) find ∇_FMR ∼0.3 in the strong correlation regime. Applying a infrared based SFR derivation method, they find somewhat lower values of ∇_FMR ∼0.2. In their study focusing on z ∼ 2.3, Sanders et al. (2018) use several example metallicity calibrations (both strong line-based/theoretical and ’direct’/empirical) and find ∇_FMR ∼0.11–0.27, hinting at potential dependence on the choice of metallicity indicator and calibration. However, as argued by Sanders et al. (2021), different metallicity calibrations (or redshift-dependent adjustments to z ∼ 0 based calibrations) may be needed to correctly derive metallicities of star-forming galaxies at different redshifts. This is taken into account in Sanders et al. (2021), where different empirical metallicity calibrations are used at z ∼ 0 and at higher redshifts. They find best-fitting ∇_FMR = 0.27 at z ∼ 0 and somewhat shallower dependence at z ∼ 2.3 (∇_FMR ∼ 0.19), although consistent within the uncertainty with z ∼ 0 determination (see their Fig. 10).

Torrey et al. (2018) show the strength of the SFR–metallicity correlation split in several log₁₀(M_*) and redshift bins as found in the IllustrisTNG simulations (see Table 1 therein). They report values in the range ∇_FMR ∼0.25 – 0.34 (except for the lowest log₁₀(M_*) ∼9 and highest mass bin log₁₀(M_*) >10.5 at z ∼ 0, where noticeably weaker ∇_FMR is found: 0.19 and 0.1, respectively), with no clear mass or redshift dependence.

We use ∇_FMR0 = 0.27 as recently found by Sanders et al. (2021) as our fiducial choice. Given the uncertainty of this parameter, we consider values between 0.17 and 0.3 dex to discuss the sensitivity of our results to this choice.⁵

3.1.6 Calculating metallicities of galaxies at different z with Z_O/H(SFR, M_*)

Z_O/H(SFR, M_*) constructed as described in this section is fully determined by the choice of ∇_FMR0 and the parameters of the z ∼ 0 MZR and SFMR. We introduce the variations of those local relations that we explore in this study in Section 3.4. We further assume that Z_O/H(SFR, M_*) does not evolve with redshift–metallicity evolution with redshift is then a result of an evolving GSMF and SFMR. We sample galaxy masses from the redshift-dependent GSMF and their SFR from the distribution centred around the redshift-dependent SFMR, as discussed in Section 2. We then use Z_O/H(SFR, M_*) to assign the metallicity. Furthermore, observational studies indicate the residual scatter around the FMR σ_FMR ∼ 0.05 dex. To account for that, we add a normally distributed scatter σ_FMR = 0.05 dex to metallicities assigned with our Z_O/H(SFR, M_*).

3.2 Differences with respect to Chruslinska & Nelemans (2019)

Even though the metallicity distribution of galaxies used in Chruslinska & Nelemans (2019) relies on the redshift-dependent MZR, the authors include the simplified form of the mass–metallicity-SFR dependence within their framework, assuming that the galaxy offsets from the average MZR and SFMR are fully anticorrelated, i.e. setting the coefficient from equation (2) to |$\nabla _{\mathrm{ FMR}}=\frac{\sigma _{\rm MZR}}{\sigma _{\rm SFMR}}\approx 0.33$|⁠, where σ_MZR = 0.1 and σ_SFMR = 0.3 dex describe the scatter around the average MZR and SFMR, respectively. The same ∇_FMR is used independent of M_*, SFR, or redshift. In the strong correlation regime discussed in Section 3.1.5 (at low/intermediate M_* and high SFR) and at z ∼ 0 their description of the mass–metallicity-SFR dependence is effectively the same as implemented in this study (except for the lower value of ∇_FMR used in this work). However, the strength of the observed SFR–metallicity correlation appears to weaken at high M_* and low SFR/sSFR. This behaviour is not captured with the simple description given by equation (2) with fixed ∇_FMR. Another important difference is that, in this study the Z_O/H(M_*,SFR) (once defined with z ∼ 0 relations) is assumed to be redshift invariant. Chruslinska & Nelemans (2019) use equation (2) to calculate offsets relative to SFMR and MZR as found at any given redshift. This means that, due to strong MZR evolution at z ≳ 2, their Z_O/H(M_*,SFR) is redshift dependent.

3.3 FMR example: comparison with Sanders et al. (2021)

In this section, we present an example Z_O/H(SFR, M_*) obtained with our phenomenological model, where we keep all the input assumptions as similar as possible to Sanders et al. (2021). We choose this study, as the authors provide estimates of all the ingredients that are necessary to construct our Z_O/H(SFR, M_*) (the z ∼ 0 MZR, SFMR, and ∇_FMR0),⁶ which allows for a direct comparison with their results. The resulting Z_O/H(SFR, M_*) is shown in Fig. 3. In the bottom panels we include two other common 2D projection of the FMR: the Z_O/H–log₁₀(SFR) plane (bottom left panel) and the Z_O/H–log₁₀(sSFR) plane (bottom middle panel). The thick grey line in the main panel of Fig. 3 shows the z ∼ 0 MZR as given in Sanders et al. (2021). The coloured points around that line represent the z ∼ 0 population of star-forming galaxies from our model.⁷ The black line indicates a fit to the maximum density region occupied by the z ∼ 0 galaxy sample described above, to ensue that the z ∼ 0 MZR is reproduced. The right bottom panel shows the Z_O/H residuals around the MZR obtained that way for several example log₁₀(M_*), which shows that the typically indicated intrinsic width of the relation σ_MZR ∼ 0.1 dex is reasonably reproduced (note that this quantity is not an input in our model). The thick solid coloured lines in the main panel show Z_O/H(SFR, M_*) for various fixed log₁₀(SFR) values. The faint coloured solid lines in the background were obtained with the best-fitting FMR given by Sanders et al. (2021) (see equation 10 therein), plotted roughly in the range of M_* and SFR probed by their galaxy sample. We also plot their z ∼ 2.3 (squares) and z ∼3.3 (circles) data points (obtained for stacked spectra of the observed galaxies, see Table 1 therein). The inner (outer) colours in all symbols correspond to upper (lower) bound on the log₁₀(SFR) within the uncertainties provided by the authors for each of the stacks. The long dashed and dotted grey lines indicate the projected z ∼ 2.3 and z ∼ 3.3 MZR, calculated in the same way as the black line, but using the z ∼ 2.3 and z ∼ 3.3 SFMR from Sanders et al. (2021), respectively, and our GSMF estimated at the corresponding redshifts. The overall agreement is remarkable.

Figure 3.

Example Z_O/H(SFR, M_*) obtained with our method (coloured lines – different log₁₀(SFR) values) and using the z ∼ 0 MZR, SFMR, and ∇_FMR0 from Sanders et al. (2021) (S20). Coloured points are sampled from the GSMF and SFMR at z ∼ 0. Black line is plotted to show that z ∼ 0 MZR (shown as thick solid grey line and shaded area – 3 σ_MZR region) is well reproduced with our Z_O/H(SFR, M_*) and z ∼ 0 GSMF and SFMR. The grey dashed lines show fits to the maximum density of galaxies sampled from GSMF and SFMR at z = 2.3 and 3.3 to indicate the projected MZR evolution. Data points are z ∼ 2.3 and 3.3 stacks from S20. Faint lines show FMR fitted by S20. Coloured dashed lines show Z_O/H(SFR, M_*) as would be obtained with equation (2) and ∇_FMR0 = 0.27 used at all masses and SFRs. Projections on to Z_O/H–log₁₀(SFR) and Z_O/H–log₁₀(sSFR) planes are shown in the bottom panels for several log₁₀(M_*) values (black lines – log₁₀(M_*) from 6 to 12). In the bottom left panel, we indicate the SFR corresponding to SFMR value for a given M_* at z = 0 (cross), 2.3 (square), and 3.3 (triangle). Red dots indicate the SFR of a galaxy ±3σ_SFMR away from the SFMR at z ∼ 0. The rightmost bottom panel shows the offsets of z ∼ 0 model galaxies from the z ∼ 0 MZR, plotted for several mass bins.

Open in new tabDownload slide

We also show Z_O/H(SFR, M_*) for various fixed log₁₀(SFR) values as would have been obtained with the approach used in Chruslinska & Nelemans (2019) i.e. using equation (2) with ∇_FMR = ∇_FMR0 = 0.27 fixed (i.e. the same at all SFRs and masses, see coloured dashed lines in the main panel of Fig. 3). It can be seen that the Z_O/H(SFR, M_*) constructed that way is identical with our model at low/intermediate M_* and high SFRs (the strong correlation regime), and starts to deviate at high M_* and low SFR (above the z ∼ 0 MZR), where ∇_FMR < ∇_FMR0. This demonstrates that the mass and SFR dependence of the SFR–metallicity correlation discussed in Section 3.1.5 is present in our description.

3.4 Model variations

In this study we consider several variations of the base z ∼ 0 MZR and SFMR relations, representing extreme choices of the shapes of the two relations and the MZR normalization Z_{O/H; MZR0} reported in the literature. We also consider two variations of the GSMF: either with fixed low mass end slope (α_GSMF = −1.45), or with α_GSMF steepening with redshift (see Section 3.1 in Chruslinska & Nelemans 2019 for the details). This allows us to explore the extremes of the f_SFR(Z,z) distribution. We briefly introduce each of the MZR and SFMR variations used to construct our Z_O/H(SFR,M_*) below.

3.4.1 MZR variations

Figure B1.

Fraction of the total SFRD that occurs below Z_O/H⊙-1 as a function of redshift/lookback time, plotted for the model variations shown in Fig. 12. Different colours correspond to different assumptions about a_MZR and starburst galaxies (see legend). Thick solid (dashed) lines mark the upper (lower) edges of the corresponding ranges for each of the model variations. Thin solid lines show the location of the upper edge obtained under the assumption of a non-evolving low-mass end of the GSMF.

Open in new tabDownload slide

Figure B2.

Fraction of the total SFRD that occurs above Z_O/H⊙ as a function of redshift/lookback time, plotted for the model variations shown in Fig. 13. Different colours correspond to different assumptions about a_MZR and starburst galaxies (see legend). Thick solid (dashed) lines mark the upper (lower) edges of the corresponding ranges for each of the model variations.

Open in new tabDownload slide

3.4.2 SFMR variations

To describe the z ∼ 0 SFMR we follow Chruslinska & Nelemans (2019) and assume the low/intermediate mass slope a_SFMR and normalization b_SFMR based on Boogaard et al. (2018). As discussed in Chruslinska & Nelemans (2019), the shape of the high-mass end is debated: while many authors report a varying degree of flattening above a certain mass (e.g. Schreiber et al. 2015; Tomczak et al. 2016; Bisigello et al. 2018), others see no evidence for the change of slope (e.g. Renzini & Peng 2015; Pearson et al. 2018). We consider the same extreme variations as described in Chruslinska & Nelemans (2019): ’no flattening’ – a linear relation between log₁₀(SFR) and log₁₀(M_*) with a single slope at all masses and ‘sharp flattening’ – with a_SFMR ∼ 0 at high masses [using the parametrization from Tomczak et al. (2016)]⁹. We focus on those two extremes to discuss our results, but use additional ‘moderate’ variations in some of the figures introducing our FMR model. Those assume that the slope of the SFMR changes from a_SFMR to a_SFMR2 < a_SFMR at log₁₀(M_*) = log|$_{10}\left(\rm \mathit{ M}_{0;SFMR0} \right)$|⁠. Variation moderate/S20 assumes a high-mass slope as in the z ∼ 0 SFMR shown in Sanders et al. (2021), while moderate/S14 follows the prescription of Speagle et al. (2014). The relevant parameters for all the SFMR variations are given in the bottom rows of Table 1. The SFMR evolution with redshift is described as in Chruslinska & Nelemans (2019).

3.4.3 The FMR for the considered model variations

In this section, we summarize the key differences between Z_O/H(SFR, M_*) obtained for the different variations of the z ∼ 0 MZR, SFMR, and ∇_FMR0 considered in this study. It can be seen from equation (6) that the impact of the SFMR parameters on the SFR-dependent turnover mass |$\rm log_{10}( \mathit{ M}_{0;SFR})$| (and on α) is reduced by the multiplication by ∇_FMR ≲ 0.3. Considering a range of a_SFMR ∼ 0.7–1 spanned by different determinations of the low-mass slope of the SFMR present in the literature (see e.g. fig. 10 in Boogaard et al. 2018) instead of using a fixed value a_SFMR = 0.83 would only affect α by ≲0.08. This is illustrated by the grey bands in the top panels in Fig. 2. At the same time, for the range of a_MZR considered in this study α varies between ∼0.2 and 0.55 (compare the orange and blue lines in Fig. 2, where smaller values correspond to steeper MZR slopes). Variation in ∇_FMR between 0.17 and 0.3 affects α by ∼0.15 (see green lines and bottom left panel in Fig. 2). The choice of MZR has a decisive role in setting the Z_O/H(SFR, M_*). ∇_FMR0 also visibly affects the relation, while SFMR has a relatively mild impact on its shape within our framework. Z_O/H(SFR, M_*) obtained for the different variations of the z ∼ 0 MZR considered in this study and for our fiducial ∇_FMR0 = 0.27 is shown in Fig. 4. The shift in normalization between the PP04/S20 and KK04a/b variations (compare top and bottom panels) as well as the difference in the spacing between the different log₁₀(SFR) lines depending on a_MZR (compare right-hand and left-hand panels) are evident. The analogous figures showing the impact of different SFMR (Fig. A3) and ∇_FMR0 (Fig. A4) choices are shown in Appendix  A. The additional Z_O/H(SFR, M_*) variation with an SFR-dependent ∇_FMR is discussed in Appendix  C and illustrated in Fig. C1. The choice of the SFMR variation (its high-mass end) has the strongest effect on the β parameter (see Fig. A2 in the Appendix  A, which shows the β – log₁₀(SFR) relation for different choices of the local scaling relations). In practice, cases with an SFMR with no flattening/only mild change of slope are well described with a single value of β. The dependence on SFR becomes apparent in cases with SFMR showing a significant deviation from a single power law. The obtained values of β are generally smaller for steeper MZR. β is only weakly affected by the choice of ∇_FMR0 (shifting towards smaller values with decreasing ∇_FMR0). We note that using the z ∼ 0 MZR from Curti et al. (2020), we can recover their best fit β ∼ 2 (when we assume single power-law SFMR). Before showing the results, we first discuss the treatment of starburst galaxies in our models.

Figure 4.

Z_O/H(SFR, M_*) obtained as described in Section 3.1 – different colours show Z_O/H(SFR, M_*) at fixed log₁₀(SFR) (increasing from left to right, see legend). Different panels correspond to different choices of the z ∼ 0 MZR (indicated by the thick grey line). All panels assume SFMR with a moderate flattening (z ∼ 0 SFMR with a_SFMR = 0.72 at high masses and a_SFMR = 0.83 at low masses) and ∇_FMR0 = 0.27. Grey dashed lines indicate the projected MZR evolution as obtained by performing fits to maximum density of galaxies in 12 + log₁₀(O/H)–log₁₀(M_*) plane, where the galaxies were sampled from GSMF (assuming α_GSMF = α_fix) and SFMR at several higher redshifts and their Z_O/H was assigned using the Z_O/H(SFR, M_*) shown in each of the panels. Horizontal orange lines indicate solar and 10 per cent solar metallicity, assuming Z_O/H⊙ = 8.83 (Grevesse & Sauval 1998).

Open in new tabDownload slide

4 STARBURST GALAXIES

A common approach to describe the SFR distribution of galaxies is to use a Gaussian distribution centred around the redshift-dependent SFMR.¹⁰ In reality, the SFR of star-forming galaxies at fixed stellar mass seems to follow a bimodal, double-Gaussian shape. The secondary peak of this distribution is attributed to starburst (SB) galaxies – strong SFMR outliers that feature SFR a few times higher¹¹ than those of the regular star-forming galaxies of the same mass and redshift. Several authors estimate the fraction of starburst galaxies (f_SB; the ratio between the number of galaxies associated with the starburst component of the SFR distribution and the total number of star-forming galaxies in the considered mass and redshift range) and report values f_SB ∼ 2–3 per cent (e.g. Rodighiero et al. 2011; Béthermin et al. 2012; Sargent et al. 2012; Ilbert et al. 2015; Schreiber et al. 2015). Despite the relatively low f_SB values, the starburst’s contribution to the total cosmic SFRD could still amount to |$\sim 10{{\ \rm per\ cent}}$| at z ∼ 2 due to their high SFR. Boco et al. (2021) use the double Gaussian distribution of galaxy SFRs from Sargent et al. (2012) and suggest that accounting for the starburst component can improve the consistency between the cosmic SFRD at z≳ 2 determined with the use of galaxy stellar mass functions paired with SFMR (as used in Chruslinska & Nelemans 2019) and that estimated with the use of SFR functions (which better account for the SFR of dusty galaxies at high redshifts; as used in Boco et al. 2019).

Crucially, the above mentioned f_SB determinations are based on galaxy samples limited to relatively massive objects (log₁₀(M_*/M_⊙) > 10) and z ≲ 2. Studies of Caputi et al. (2017) and Bisigello et al. (2018) extend the analysis of the distribution of star-forming galaxies in the log₁₀(SFR)–log₁₀(M_*) plane to much lower masses and higher redshifts. While their results are consistent with the previous determinations at high log₁₀(M_*), they show that f_SB is a strong function of stellar mass (increasing towards lower log₁₀(M_*)) and z (increasing with redshift). Moreover, they indicate that starburst galaxies follow a distinct sequence in the log₁₀(SFR)–log₁₀(M_*) plane that is located ∼1 dex above the SFMR. This offset of the starburst sequence relative to SFMR is considerably higher than previously reported (e.g. Béthermin et al. 2012; Sargent et al. 2012). This suggests that the contribution of starburst galaxies to the total SFRD budget can be much higher than previously estimated. If starburst galaxies follow the general FMR (as suggested by the results of Hunt et al. 2012; to our knowledge, there is no evidence to the contrary), they would contribute to the star formation at relatively low metallicities compared to galaxies on the SFMR. Therefore, they affect the low metallicity tail of the f_SFR(Z,z) – crucial for the discussion of the origin of transients as long gamma-ray bursts and double black hole mergers. We aim to discuss the possible impact of starbursts on the f_SFR(Z,z) in view of the results reported in Caputi et al. (2017) and Bisigello et al. (2018). In the following Section 4.1, we describe the method used to include the contribution of starburst galaxies within our framework.

4.1 Method and considered variations

To account for the contribution of starburst galaxies in our calculations, we follow the procedure outlined below:

• At each redshift, we sample log₁₀(M_*) of star-forming galaxies from the galaxy stellar mass function as described in Chruslinska & Nelemans (2019).

• We use f_SB to describe the fractions of starburst galaxies and regular star-forming galaxies at each log₁₀(M_*) and z. Our choice is outlined further in this section.

• The SFR of regular galaxies is given by the SFMR from Chruslinska & Nelemans (2019). The SFR of starburst galaxies at each log₁₀(M_*) follows a normal distribution with scatter σ_SB. The peak of the distribution can be related to the galaxy’s log₁₀(M_*) with a linear relation: |$\rm log_{10}(SFR_{SB}) = a_{SB} log_{10}( \mathit{ M}_{*}) + b_{SB}$|⁠, to which we refer as the starburst sequence We set the value of b_SB by defining the offset of the starburst sequence from the SFMR, i.e. b_SB = b_SFMR + Δ_{SB − SFMR}. The assumed parameters are given further in this section.

• To describe the metallicity at which the starburst galaxies produce stars, we assume that they follow the same FMR as regular star-forming galaxies (described in Section 3.1).

We consider two sets of parameters describing the properties of starbursts. First, we follow the same implementation as used in the recent study by Boco et al. (2021), which is based on the works of Sargent et al. (2012) and Béthermin et al. (2012). This implementation assumes a constant f_SB = 0.03 (independent of mass and redshift), a starburst sequence with scatter σ_SB = 0.24, located Δ_{SB − SFMR} = 0.59 dex above the SFMR and parallel to the SFMR (specifically, we assume the starburst sequence slope of a_SB = 0.83, corresponding to the low-mass slope of the SFMR from Chruslinska & Nelemans 2019,¹²). In light of the results of Caputi et al. (2017) and Bisigello et al. (2018) this implementation severely underestimates both f_SB and the SFR of galaxies on the starburst sequence. Therefore, it likely provides the absolute lower limit on the contribution of starbursts. Secondly, we follow the results of Caputi et al. (2017) and Bisigello et al. (2018) to model the mass and redshift dependence of f_SB and the properties of the starburst sequence. We provide the details of this implementation in Sections 4.1.1 and 4.1.2 and refer to it as the B18/C17 implementation in the reminder of this paper.

4.1.1 The fraction of starbursts

Figure 5.

Fraction of starbursts as a function of galaxy stellar mass. Thick lines show the assumed dependence at z ∼ 0.75, 1.5, and 2.5 based on the data from Bisigello et al. (2018) and at z ∼ 4.4 based on Caputi et al. (2017). Thin lines show the interpolation between those redshifts and extrapolation to z = 0 (thin horizontal dashed line), plotted with 0.25 step in redshift at z < 1.5 and with 0.5 step in redshift above. Vertical dashed lines mark log₁₀(M_fix) (see Table 2) for different redshift bins. Open symbols mark estimates where the galaxy sample from Bisigello et al. (2018) is below 90 per cent mass completeness.

Open in new tabDownload slide

All stellar masses in equation (8) are in solar units M_⊙. The adopted coefficients are given in Table 2. We assume that the above relation holds strictly in the middle of each redshift bin and interpolate between them to describe f_SB at 0.5 < z < 4.4. At z > 4.4 we use the relation from z= 4.4. That way we obtain a conservative estimate of the starburst contribution at higher redshifts. We extrapolate the constructed dependence to lower redshifts, setting a constant starburst fraction of 3 per cent at all masses at z = 0. At z = 0.75, 1.5, and 2.5, log₁₀(M_fix) is the lower edge of the lowest mass bin where the galaxy sample from Bisigello et al. (2018) is within 90 per cent stellar-mass completeness. We make a conservative assumption and use a fixed f_SB value below that mass (i.e. |$f_{\mathrm{ SB}\mathrm{ }} =\rm a_{1} log_{10}( \mathit{ M}_{\rm fix}) + b_{1}$|⁠). At z = 4.4, log₁₀(M_fix) = 9.25 corresponds to the lowest mass for which the fraction of starbursts has been estimated in Caputi et al. (2017). At this redshift, f_SB is fixed at the same mass as at z = 2.5 (i.e. log₁₀(M_*) = 8.24) and we assume a linear relation to describe the dependence between log₁₀(M_*) = 8.24 and log₁₀(M_*) = 9.25. This exception is necessary to avoid a decrease in the fraction of starbursts at log₁₀(M_*) < 9.25 at z > 2.5, which would break the trend seen in the data (see Fig. 5). At log₁₀(M_*) > 10.5 there is a lot of scatter in the data and the trend is inconclusive. In this mass range we use a fixed value f_SB = 0.03 at all redshifts, as found in earlier studies that focused on the most massive galaxies (e.g. Rodighiero et al. 2011; Schreiber et al. 2015). The steeper, lower mass part of f_SB–log₁₀(M_*) dependence is well described with the same slope across all four redshift bins. Therefore, for simplicity we assume a₁ = −0.3. The slope of the relation in the higher mass part (a₂) and the mass log₁₀(M_{0; SB}) separating the high and low mass parts of the relation increase with redshift.

4.1.2 The starburst sequence

We guide our description of the SFR distribution of starburst galaxies with the results shown in figs 3 and 7 from Caputi et al. (2017) and figs 6 and 7 from Bisigello et al. (2018). The resulting starburst sequence is shown in Fig. 6. The width of the starburst sequence in Caputi et al. (2017) and in Bisigello et al. (2018) is smaller that that of the SFMR, although the values are not given. Guided by the results shown in the figures, we assume σ_SB = 0.2 dex. Both Caputi et al. (2017) and Bisigello et al. (2018) find a starburst sequence that is steeper than the SFMR. There is no clear evidence for a_SB evolution with redshift. For simplicity, we assume a_SB = 0.94 (average between the best-fitting values in three redshift bins from Bisigello et al. and high redshift estimate from Caputi et al.). To set the offset of the starburst sequence from the SFMR Δ_{SB − SFMR} we focus on the results for the intermediate masses log₁₀(M₁₀/M_⊙) ∼ 9-9.5 (where the sample is complete and the SFMR is not affected strongly by the potential flattening at the high-mass end). We find that Δ_{SB − SFMR} in both Bisigello et al. (2018) and Caputi et al. (2017) is about 1 dex and we assume this value in our calculations. Note that the SFMR shown in Caputi et al. (2017) is likely an upper limit on the SFMR location at z ∼ 4.4, which suggests that at those high redshifts the offset might be even larger.

Figure 6.

Starburst sequence (SB – dashed lines) compared with the SFMR for extreme assumptions about the high-mass slope of the SFMR: no flattening (thick solid lines) and sharp flattening (thin solid lines). Colours correspond to different redshifts. Points sampled from the GSMF and assigned to the SFMR or starburst sequence using starburst fraction shown in Fig. 5 are plotted for z = 0 and z = 5 to illustrate the increasing fraction of star-forming galaxies that populate the starburst sequence at higher redshifts and lower stellar masses. The properties of starbursts are based on the results from Caputi et al. (2017) and Bisigello et al. (2018).

Open in new tabDownload slide

5 RESULTS: THE DISTRIBUTION OF THE COSMIC SFRD OVER METALLICITIES AND REDSHIFT

In this section, we discuss the distributions of the cosmic SFRD over metallicities and redshift (f_SFR(Z,z)) obtained for different variations of our observation-based model (including different choices of the local MZR, SFMR, GSMF, ∇_FMR0, and prescriptions to account for starburst galaxies). Ultimately, we aim to explore the extreme f_SFR(Z,z) cases in terms of the amount of SFRD occurring at low (below 10 per cent solar metallicity; Z_O/H ≤ Z_O/H⊙−1) and high (above solar Z_O/H ≥ Z_O/H⊙) metallicity. In the remainder of this paper we distinguish between two sets of models based on the choice of the local scaling relations (MZR and SFMR): the ‘low metallicity’ f_SFR(Z,z) cases – obtained with z ∼ 0 MZR with low normalization (PP04 or S20) and SFMR with sharp flattening at high masses, and the ‘high metallicity’ f_SFR(Z,z) cases – obtained with z ∼ 0 MZR with high normalization (KK04a or KK04b) and SFMR with no flattening. Other combinations of the local MZR and SFMR lead to more moderate metallicity distributions.

We note that in order to describe the f_SFR(Z,z) at high redshifts and low metallicities (as shown in this Section and needed, for instance, in applications to gravitational wave astrophysics) one needs to extrapolate the FMR well beyond the regions where it is constrained by current observations (in particular to z ≳ 3 and M_* ≲ 10⁸ M_⊙). The importance of the assumed z ≳ 3 extrapolation is discussed in Section 5.1, where we compare the f_SFR(Z,z) obtained with the non-evolving FMR and with the redshift-dependent MZR based approach used in Chruslinska & Nelemans (2019). Metallicity of M_* ≲ 10⁸ M_⊙ galaxies assigned with our FMR is primarily sensitive to the MZR slope and normalization and the strength of the SFR–metallicity anticorrelation at low masses/SFRs. Those factors are discussed in Section 5.1.1 and further in Appendix  C. In Section 5.2, we discuss the f_SFR(Z,z) under different assumptions about the contribution of starburst galaxies.

5.1 f_SFR(Z,z) with redshift-invariant FMR

The comparison between the f_SFR(Z,z) constructed with the use of redshift-invariant Z_O/H(SFR,M_*) and the corresponding distributions from Chruslinska & Nelemans (2019) are shown in Fig. 7. We compare the variations with z ∼ 0 MZR and SFMR choices as in the high (left) and low (right) metallicity extremes from Chruslinska & Nelemans (2019) (see Section 4.2 and Fig. 8 therein). In this section, we do not explicitly include starburst galaxies, so that all assumptions that are not related to the metallicity distribution of galaxies are the same as in Chruslinska & Nelemans (2019).

Figure 7.

Distribution of the star formation rate density (SFRD) at different metallicities and redshift (z) in this paper (background colours, brown contours; assuming non-evolving FMR and ∇_FMR0 = 0.27) compared to Chruslinska & Nelemans (2019) (white contours, assuming redshift-dependent MZR). Colour indicates the amount of SFRD happening in different redshift and metallicity bins, contours indicate constant SFRD and are plotted for 0.01, 0.05, and 0.1 [M_⊙/Mpc³yr] (with the highest value corresponding to the most inward contour). Left – model variation with KK04a z ∼ 0 MZR and SFMR with no flattening at high masses compared with the high metallicity extreme from Chruslinska & Nelemans (2019). Right – model variation with PP04 z ∼ 0 MZR and SFMR with sharp flattening at high masses compared with the low metallicity extreme from Chruslinska & Nelemans (2019). All models assume GSMF with non-evolving low-mass end. Orange horizontal dashed lines indicate solar, 10 per cent solar and 1 per cent solar metallicity (assuming Grevesse & Sauval 1998 solar metallicity scale). The right metallicity axis was obtained from Z_O/H assuming solar abundance ratios. Note that beyond z ≳ 3 the distribution relies on extrapolation of the MZR/FMR evolution with redshift – in this respect, the plotted contours contrast somewhat extreme assumptions.

Open in new tabDownload slide

Figure 8.

Example illustrating the effect of different assumptions about the low mass end slope of the GSMF (α_GSMF) on the f_SFR(Z,z). Top panel: brown constant SFRD contours and background colours – variation with α_GSMF steepening with redshift (i.e. predicting increasing number density of low-mass galaxies), white contours – fixed α_GSMF = −1.45. Contours are plotted for 0.01, 0.05, and 0.1 [M_⊙/Mpc³yr]. Both model variations assume ∇_FMR0 = 0.27 and the same SFMR and z ∼ 0 MZR (indicated in the figure). f_SFR(Z,z) significantly deviate at z ≳ 3: a redshift-dependent α_GSMF leads to most of the SFRD happening at low metallicity and to a higher total SFRD(z). Bottom panel: ratio of the total SFRD at each redshift as obtained in the two variations.

Open in new tabDownload slide

Fig. 7 shows that the metallicity distributions start to deviate around z ≳ 1.5 (compare brown and white contours), where the redshift-dependent MZR predicts a steeper decrease in metallicity than what results from the non-evolving FMR. The difference becomes striking at z ≳ 3, but we stress that neither the MZR nor the FMR is currently constrained in this redshift regime. In particular, there is no guarantee that the FMR holds or continues to show weak/no redshift evolution beyond that redshift. On the other hand, as discussed in Section 3, the extrapolated MZR evolution as assumed by Chruslinska & Nelemans (2019) is likely to overestimate the rate of decrease in metallicity. Therefore, in Fig. 7 we contrast two seemingly extreme assumptions. Until better observational constraints are available, this comparison can serve to illustrate (likely a conservative) range of uncertainty of the high redshift part of f_SFR(Z,z) resulting from the extrapolated evolution of the galaxy metallicity distribution with redshift. The extrapolated MZR evolution leads to a z ∼ 10 peak metallicity that is almost 2 dex lower than what results from the non-evolving FMR assumption. The latter leads to SFRD concentrated at higher metallicities (irrespective of the model variation), but the extended low metallicity tail is still present at all redshifts. We note that the difference between the f_SFR(Z,z) obtained with a redshift dependent MZR and with the non-evolving FMR was recently discussed in Boco et al. (2021) – qualitatively our results are the same as discussed therein. However, rather than discussing a example f_SFR(Z,z) obtained with a particular FMR or MZR taken from the literature, here we model the FMR consistently with the choice of the local SFMR and MZR and can explore the uncertainties of the final result. All variations shown in Fig. 7 assume a GSMF with non-evolving low-mass slope. As discussed in Chruslinska & Nelemans (2019) (see Section 4.1 therein), this assumption has relatively little effect on f_SFR(Z,z) at z ≲ 3, but strongly affects the result at higher redshifts (both its low metallicity tail and the total SFRD). This is illustrated in Fig. 8 – the effect of the steepening low-mass end of the GSMF on f_SFR(Z,z) is analogous for other model variations.

5.1.1 The effect of the z ∼ 0 MZR slope and ∇_FMR0

The effect of the choice of the z ∼ 0 MZR slope (a_MZR) and ∇_FMR0 on f_SFR(Z,z) is shown in Fig. 9. Similarly to Fig. 7, the left-hand and right-hand panels show high and low metallicity model variations, respectively. Top (bottom) panels show variations with a_MZR ∼0.3 (a_MZR ∼0.6). Different contours illustrate the impact of ∇_FMR0.

Figure 9.

Distribution of the star formation rate density (SFRD) at different metallicities and redshift (z) for different assumptions about the z ∼ 0 MZR slope (top panels: variations with a_MZR ∼ 0.3, bottom panels: a_MZR ∼ 0.6) and ∇_FMR0 (compare different contours; background colours are plotted for variations with ∇_FMR0 = 0.27). Left-hand (right) panels show model variations with z ∼ 0 MZR and SFMR that maximize the SFRD at high (low) metallicity. Constant SFRD contours are plotted for 0.005, 0.01, 0.05, and 0.1 [M_⊙/Mpc³ yr] (the highest value corresponds to the innermost contour). All model variations assume a GSMF with non-evolving low-mass end. Variations with flatter z ∼ 0 MZR lead to a more compact metallicity distribution, but their peak decreases with redshift faster than in variations with a steeper z ∼ 0 MZR.

Open in new tabDownload slide

As can be expected, f_SFR(Z,z) for variations with steeper z ∼ 0 MZR (higher a_MZR values) have a more extended low metallicity tail at each redshift. However, variations with lower a_MZR values feature a much steeper metallicity evolution at z ≲ 4. At higher redshifts the metallicity peak of f_SFR(Z,z) decreases at a similar rate for the variations with the same ∇_FMR0 and same total SFRD (i.e. the same SFMR and GSMF assumptions), irrespective of a_MZR.¹³ This can be understood by looking at the comparison of Z_O/H(SFR,M_*) obtained for different z ∼ 0 MZR, shown in Fig. 4. In this plane, galaxies with the same log₁₀(M_*) at higher redshifts shift to higher log₁₀(SFR) lines. The rate of this shifting is dictated by the redshift-dependent SFMR. The offset in Z_O/H between the different log₁₀(SFR) curves at fixed log₁₀(M_*) is bigger for variations with lower a_MZR (compare left-hand and right-hand panels in Fig. 4) and similarly for higher ∇_FMR0 (see Fig. A4 in the Appendix  A). This translates into steeper decrease of the f_SFR(Z,z) metallicity peak with redshift in the variations with lower a_MZR and higher ∇_FMR0. At sufficiently high redshifts (where almost all galaxies occupy the high SFR, linear regime of Z_O/H(SFR,M_*)) the spacing in Z_O/H for different log₁₀(SFR) curves is set by the choice of ∇_FMR0. Lower ∇_FMR0 result in weaker metallicity evolution (compare white and brown contours in Fig. 9). Therefore, for fixed ∇_FMR0 and the same assumptions about the SFMR, the rate of metallicity evolution at z ≳ 4 and fixed log₁₀(M_*) is essentially the same. However, the full f_SFR(Z,z) distribution is also affected by the GSMF. If the low-mass end of the GSMF is allowed to steepen with redshift, the contribution of the low mass (and so – low metallicity) galaxies to the total SFRD at high redshifts is considerably higher than if the low-mass end of the GSFM is fixed (see example shown in Fig. 8). The FMR is not constrained for low-mass, low-SFR galaxies. Results discussed above assume that those galaxies are well described by the same relation as the galaxies with higher M_* and SFR, and that the strength of the SFR–metallicity correlation does not change in this part of the parameter space (i.e. ∇_FMR = ∇_FMR0). In Appendix  C, we additionally discuss the variation with an SFR-dependent ∇_FMR in which the SFR–metallicity correlation is assumed to disappear (i.e. the FMR breaks down) at SFR corresponding to z ∼ 0 SFMR galaxies with M_* ≲ 10⁸ M_⊙. Overall, this assumption has minor impact on the estimated f_SFR(Z,z) distribution (see Fig. C2).

5.2 f_SFR(Z,z): the impact of starbursts

The comparison of f_SFR(Z,z) obtained with and without including starbursts is shown Fig. 10. Top and middle panels show the same z ∼ 0 MZR and SFMR variations as in Fig. 9. Different contours correspond to different starburst implementations (introduced in Section 4.1). The effect of including starbursts in our calculations is two-fold: first, a fraction f_SB of star-forming galaxies is now assigned a higher SFR than what would result from the SFMR. This means that the total SFRD at each redshift is higher than in the model that does not explicitly account for starbursts (see bottom panel in Fig. 10). Secondly, assuming that starbursts follow the FMR, they contribute to star formation happening at relatively low (for their mass and redshift) metallicity. Therefore, including starbursts shift the peak of the f_SFR(Z,z) to lower metallicities and broadens the low metallicity part of the distribution.

Figure 10.

Distribution of the SFRD at different metallicities and redshift (z) for different assumptions about starburst galaxies: no starbursts (background colours, brown contours), starburst sequence, and f_SB model based on the results of Bisigello et al. (2018) and Caputi et al. (2017) (B18/C17, white contours) and as in Boco et al. (2021) (black dashed contours). Top and middle panels correspond to different z ∼ 0 MZR and SFMR variations (same as in Fig. 9). Contours indicate constant SFRD and are plotted for 0.01, 0.05, and 0.1 [M_⊙/Mpc³yr] (with the highest value corresponding to the innermost contour). Note that different assumptions about starburst galaxies lead to different total SFRD at each redshift – bottom panels show the ratio of the SFRD integrated over metallicities at each z in the model with starbursts (brown – B18/C17 implementation, black – Boco et al 2021 implementation) to the corresponding model without starbursts (the choice of MZR does not affect the total SFRD). All models assume GSMF with non-evolving low-mass end and ∇_FMR0 = 0.27.

Open in new tabDownload slide

Both effects are clearly seen when the variations including the B18/C17 starburst implementation (white contours in Fig. 10) are compared with the corresponding variations that do not include starbursts (brown contours). The lower edges of white contours extend to lower metallicities with respect to no starburst case at all redshifts, while the upper, high metallicity edges remain the same. This difference increases with redshift due to increasing f_SB. The total SFRD is about 2.5 times higher at z ≳ 4 in the case with starbursts (brown line, bottom panel in Fig. 10). Note that we fix f_SB beyond z = 4.4 – the highest redshift bin covered by the data in Caputi et al. (2017). If the trend seen at lower redshift continues, the difference at high redshifts would be even higher. If instead we follow the starburst implementation as used in the recent study by Boco et al. (2021) (black dashed contours), the difference with respect to cases without starbursts is negligible. This is expected, given the low fixed f_SB = 3 per cent and the fact that in this prescription the starburst sequence is within the 2σ_SFMR scatter of the SFMR. The difference is slightly more pronounced if SFMR with a sharp flattening at high masses is used (right-hand panels in Fig. 10), as in those cases there is a larger difference in SFR between massive galaxies on the SFMR and on the starburst sequence. The inclusion of starbursts as in Boco et al. (2021) also barely affects the total SFRD (see black dashed line in the bottom panel in Fig. 10).

6 METALLICITY-DEPENDENT COSMIC SFH

In this section we discuss the cosmic SFH – SFRD integrated over all metallicities as a function of redshift/cosmic time, as well as the high and low metallicity cuts of the cosmic SFH (i.e. SFRD occurring above and below certain metallicity thresholds) obtained for different model variations. Note that this is essentially a different way to present our results, which provides less detailed information than the full f_SFR(Z,z) distributions shown in Section 5. It allows us to zoom into the interesting parts of the distribution and demonstrate the uncertainty of its various parts more clearly. However, we stress that the discussed cuts of the cosmic SFH are sensitive to the considered low/high metallicity thresholds. In some cases, factors that have minor impact on the overall f_SFR(Z,z) distribution may lead to considerable uncertainty in the SFH cut for a particular choice of metallicity threshold. The 10 per cent solar/solar Z_O/H thresholds used here were chosen to zoom into the low/high metallicity tails of the distribution. In general, the relevant thresholds may vary depending on the considered problem. When showing the cosmic SFH and its various metallicity cuts, in each case we indicate the range of possible outcomes spanned by the model variations with extreme assumptions about the considered factors. We indicate which of the considered factors drive the uncertainty of the f_SFR(Z,z) in different regimes (i.e. low/high metallicity, low/high redshift). As discussed in Section 5, the difference between the f_SFR(Z,z) obtained for variations with starburst implementation as in Boco et al. (2021) (i.e. with fixed, small fraction of starbursts) and the corresponding variations with no starbursts are negligible. Therefore, in this section we only discuss the former.

6.1 Cosmic SFH

The cosmic SFH obtained for different model variations considered in this study is shown in Fig. 11. Coloured ranges distinguish different assumptions about the starburst galaxies. They span between the extreme cases resulting from different assumptions about the high-mass end of the SFMR and the low-mass end of the GSMF. Note that the assumptions about the metallicity (z ∼ 0 MZR, redshift evolution of MZR/non-evolving FMR, ∇_FMR) have no impact on the (total) cosmic SFH, but affect its different metallicity cuts. For the convenience of use and discussion, we approximate the results shown in Fig. 11 with a broken power law: SFRD = A (1 + z)^κ, fitting A and κ in several redshift ranges for the upper and lower edge of each variation. The resulting coefficients are given in Table 3. Coefficients obtained when only variations with non-evolving low-mass end of the GSMF are considered are provided in Table B1 in the Appendix  B. Considering variations with fixed assumptions about the starburst galaxies, we note the following:

• The uncertainty in the cosmic SFH at z < 2 is dominated by the assumptions about high-mass end of the SFMR: the lower edges of all ranges correspond to variations with SFMR with sharp flattening at high masses (as in low metallicity variations), while the upper edges correspond to variations with SFMR with no flattening (as in high metallicity variations). Assumptions about the GSMF have secondary role in this redshift range.

• The evolution of the low-mass end of the GSMF dominates the uncertainty at higher redshifts: the lower (upper) edges of all ranges correspond to variations with non-evolving (steepening) low-mass end of the GSMF. At the same time, the importance of the high-mass end of the SFMR is reduced due to decreasing number density of the most massive galaxies (see e.g. fig. 3 in Chruslinska & Nelemans 2019).

• The upper edges of all ranges feature an upturn in the SFRD evolution at z ≳ 4, before the sharp decrease at z > 7. This is due to strongly increasing number density of low-mass galaxies in the variations in which the low-mass end of the GSMF steepens with redshift. We note that such upturn only appears when the SFMR and GSMF are extrapolated to log₁₀(M_*) < 8 (see discussion in Chruslinska & Nelemans 2019).

• All variations show a sharp decrease in the SFRD at z ∼ 7. This is due to the rapid evolution of the GSMF normalization between z ∼ 7 and 8, seen when observational GSMF estimates from different redshifts are combined (see fig. 3 in Chruslinska & Nelemans 2019).

Figure 11.

Star formation rate density as a function of lookback time/redshift. The plotted ranges span between the estimates obtained for different variations of model assumptions about the SFMR and GSMF. The ranges in different colours correspond to different assumptions about the starburst galaxies (turquoise: as in Boco et al. (2021), brown: based on Caputi et al. (2017) – C17 and Bisigello et al. (2018) – B18). The orange region shows the corresponding range from Chruslinska & Nelemans (2019) (with no starburst contribution). The thick solid orange and brown lines show how the upper edge of the corresponding range would look like under the assumption of a non-evolving low-mass end of the GSMF (the effect on the turquoise range follows closely the one that can be seen for Chruslinska & Nelemans 2019, the lower edge is not affected by this assumption). The dashed grey line shows the estimate from Madau & Dickinson (2014) for reference. We overplot the z ≲ 2 infrared-based estimates complied by Casey et al. (2018) and the recent observational high redshift SFRD estimates by Loiacono et al. (2021), Khusanova et al. (2021), and Gruppioni et al. (2020) that account for the contribution of heavily dust obscured galaxies. Note that the latter two provide a lower limit. See text for the discussion of the comparison with observations. Where necessary, we correct for the difference in SFR due to varying IMF assumptions.

Open in new tabDownload slide

The first three effects are best seen by comparing the thick solid lines in Fig. 11 (corresponding to variations with SFMR with no flattening and GSMF with non-evolving low-mass end) with the upper edges of the relevant ranges. Comparing the variations including starbursts (turquoise/brown ranges) with those from Chruslinska & Nelemans (2019) (that do no explicitly account for their contribution; orange range), one can see that:

• The SFRD is increased in the variations including starbursts; starbursts implementation used in Boco et al. (2021) (turquoise) leads to negligible difference, while the B18/C17 implementation (brown) leads to a factor of ∼2.5 increase at high redshifts (see Section 5.2).

• Variations including starbursts span a narrower range in cosmic SFH (the lower edge is lifted with respect to no starbursts case). This is due to the fact that assigning high SFR (compared to SFMR values) to a fraction of galaxies (starbursts) effectively reduces the difference between the SFRD in model variations with no/sharp SFMR flattening. The effect is stronger at high redshifts for the B18/C17 prescription, as f_SB increases with redshift.

• the increasing fraction of starbursts in the B18/C17 variations leads to a broader peak of the cosmic SFH and a shallower SFRD decrease at high redshifts.

In Fig. 11, we contrast the total SFRD that results from various realizations of our observation-based model with several observational determinations of this quantity. At z < 2, there is about a factor of 2 (2.5) offset between the SFRD obtained with model variations that assume SFMR with no flattening at high masses (and C17/B18 starburst implementation) and the commonly used, best-fitting estimate from Madau & Dickinson (2014) (grey-dashed line). However, we note that this offset is within the scatter and uncertainty of the z ≲ 2 observational estimates – in particular, infrared based SFRD determinations suggest higher SFRD than the estimate obtained by Madau & Dickinson (2014), more in agreement with our SFRD determination based on the GSMF, SFMR, and starburst sequence (see also fig. 1 in Casey et al. 2018 – the infrared-based data complied by those authors are also included in Fig. 11; and fig. 2 and discussion in Boco et al. 2021). At higher redshifts (2 < z < 7), our model variations (even those not including starbursts) in general show shallower SFRD decrease than the estimate provided by Madau & Dickinson (2014) (SFRD ∝ (1 + z)^−2.9). This is striking for the upper edges of our estimates (see also Table 3). Shallower SFRD evolution than estimated by Madau & Dickinson (2014) is in line with the most recent high redshift observational estimates (e.g. Gruppioni et al. 2020; Khusanova et al. 2021; Loiacono et al. 2021, although the uncertainties in those measurements are still substantial), indicating that a significant fraction of the SFRD at z ≳ 2 still occurs in dust obscured galaxies (and hence was missed by the earlier – mostly UV light-based surveys). Those estimates¹⁴ are also included in Fig. 11. However, the comparison of our results with those estimates is not straightforward for two main reasons: (a) observations may be incomplete: the estimate provided by Khusanova et al. (2021) is based on a UV selected sample and may not fully account for extremely dusty galaxies, while the estimate provided by Gruppioni et al. (2020) does not account for the contribution of UV sources (b) the faint (low mass) end of the galaxy population is not directly probed by observations. All estimates shown in Fig. 11 include (different) extrapolations to account for their contribution to the total SFRD. In our models we extrapolate all of the empirical relations down to log₁₀(M_*/M_⊙) = 6 – as discussed above, the result of this extrapolation at high z is very sensitive to the GSMF slope. Extrapolations used in the high redshift observational estimates give less weight to low-mass galaxies than our models with steepening low-mass end of the GSMF. The estimate by Khusanova et al. (2021) includes extrapolations down to M_*/luminosity limit comparable to the extrapolation limit used in our study, but is averaged over the results obtained with different GSMF/UV luminosity functions (with different low mass end slopes) from the literature. Gruppioni et al. (2020) extrapolate the infrared luminosity function down to low luminosities. However, such galaxies are expected to be relatively unobscured and the faint-end slope of the infrared luminosity function is much flatter than that of the UV luminosity function (and the GSMF used in our study) at high redshifts – in that sense, this estimate underestimates the contribution of low-mass galaxies to the total SFRD. Therefore, the observational estimates are likely lower limits when compared to our models. We note that all model variations with fixed low mass end of the GSMF and no/small fraction of starbursts fall below the recent z ∼ 3 limit from Gruppioni et al. (2020). Given the fact that this determination underestimates the contribution of low-mass galaxies to the total SFRD with respect to our models, this offset may hint at the need for higher contribution of massive galaxies than included in those model variations (e.g. higher fraction of starburst galaxies or higher overall normalization of the SFMR at those redshifts). We also note that the z ∼ 4.5 estimate by Khusanova et al. (2021) falls above our estimates obtained with fixed low mass end of the GSMF and no/small fraction of starbursts, but it is unclear how to properly compare this averaged estimate with our results. All model variations are consistent with the limits obtained by Gruppioni et al. (2020) and Loiacono et al. (2021) (although the latter is not particularly constraining due to large uncertainty that overlaps with all of the considered model variations) at similar redshifts.

6.2 A low metallicity cut of the cosmic SFH

6.2.1 Note on the low metallicity threshold

When referring to low metallicity, throughout this paper we use a limit of <10 per cent solar oxygen abundance (with Z_O/H⊙ following Grevesse & Sauval 1998) and the numbers quoted later in this section necessarily depend on this choice. This threshold is chosen arbitrarily and serves for illustrative purposes. In general, the relevant definition of low metallicity varies depending on the discussed problem and application. For instance, the abrupt drop commonly seen in the metallicity dependence of the formation efficiency of black hole – black hole mergers formed in the course of isolated binary evolution appears at metallicities 10–50 per cent solar (depending on the considered model of binary/stellar evolution, e.g. Chruslinska et al. 2019; Broekgaarden et al. 2021). In that sense, depending on the model, the interesting low metallicity threshold could be as high as 50 per cent solar metallicity. Furthermore, from the point of view of evolution of massive stars and the related transients, the most relevant metallicity tracer is iron rather than oxygen (which serves as a metallicity proxy in this study; the oxygen abundance is the most easily obtained from emission line measurements and so it is the most commonly used metallicity proxy in studies of interstellar medium of galaxies). That is primarily due to high sensitivity of the radiation-driven mass-loss rates of hot, massive stars to the abundance of this element (e.g. Vink, de Koter & Lamers 2001; Vink & de Koter 2005; Gräfener & Hamann 2008; Sanders et al. 2021). In general, interstellar medium (ISM) of galaxies is enriched in oxygen and iron on different time-scales due to the fact that different sources dominate the production of those elements (core-collapse supernovae with ≲10 Myr delay with respect to star formation for oxygen and SNIa supernovae with ∼0.1–1 Gyr delay for iron; e.g. Wheeler, Sneden & Truran 1989). This means that especially young, highly star-forming galaxies (i.e. typical at high redshifts) are expected to show overabundance of oxygen relative to iron with respect to solar ratios. Recent studies indeed find evidence for supersolar O/Fe≳2 (O/Fe)|_⊙ in the ISM of z ∼ 2 star-forming galaxies (e.g. Steidel et al. 2016; Topping et al. 2020; Cullen et al. 2021). Consequently, oxygen (or, more generally, α-element)-enhanced abundance ratios are also found in local stars that formed at earlier cosmic epochs (e.g. Bensby, Feltzing & Lundström 2004; Tolstoy, Hill & Tosi 2009). There is no universal way to translate metallicity measure based on oxygen abundance to iron abundance, in particular the typical oxygen to iron ratio in the ISM (and the interesting low metallicity threshold – if oxygen is used as a proxy) can be expected to vary with redshift and possibly galaxy characteristics. A more thorough discussion of this issue is beyond the scope of this work and is left to future studies.

6.2.2 Results

SFRD that occurs at Z_O/H < Z_O/H⊙-1 as a function of redshift/time for different variations of our model is shown in Fig. 12. Note that the assumptions about the high-mass end of the SFMR have no impact on the low metallicity tail of f_SFR(Z,z) and the low metallicity SFH. Similarly as for the cosmic SFH, we provide the power-law fits to the upper and lower edges of the ranges shown in Fig. 12 fitted in several z ranges. Those are listed in Tables B3 and B2 in the Appendix  B. The following features are common to all model variations:

• At z ≲ 2, the SFRD occurring at low metallicity increases faster than the total SFRD. The exact slope depends on the variation. In this redshift range the width of the ranges is set by the assumptions about the z ∼ 0 MZR: lower/upper edges of all ranges correspond to variations with high/low MZR normalization (as in the high/low metallicity variations). At those redshifts, assumptions about the low-mass end of the GSMF have negligible impact on the amount of star formation happening at low metallicity.

• Beyond the peak of the cosmic SFH assumptions about the low-mass end of the GSMF set the width of the ranges and the behaviour of the upper edges. The lower/upper edges correspond to variations with high/low MZR normalization and non-evolving/steepening GSMF low mass end.

• If α_GSMF is fixed, the low metallicity cut of the SFH shows a broad peak between z = 2–4 (depending on the variation) and a relatively mild evolution up to z ∼ 7.

• If α_GSMF evolves with redshift, the increasing abundance of low-mass galaxies (forming stars at low metallicity) causes a continuous increase in the low metallicity cut of the SFH up to z ∼ 7.5–8.

• All variations show a sharp decrease in the low metallicity SFRD at z ≳ 8, reflecting the sharp drop in the total SFRD at those redshifts.

Figure 12.

Star formation rate density at low metallicity [i.e. below 10 per cent solar, assuming solar oxygen abundance Z_O/H⊙ = 8.83 Grevesse & Sauval (1998)] as a function of lookback time/redshift. The plotted ranges show estimates obtained for different variations of assumptions about the starburst (SB) galaxies (blue/turquoise: as in Boco et al. (2021), brown/beige: based on Caputi et al. (2017) – C17 and Bisigello et al. (2018) – B18) and low-mass slope of the redshift z ∼ 0 MZR (a_MZR). The orange region shows the range of the low metallicity SFHs obtained by Chruslinska & Nelemans (2019). Thick solid lines in different colours indicate the location of the upper edge of the corresponding ranges obtained with the assumption of non-evolving low-mass end of the GSMF, demonstrating the importance of this assumption beyond the peak of the cosmic SFH. All variations from this study assume ∇_FMR0 = 0.27.

Open in new tabDownload slide

The impact of the assumed z ∼ 0 MZR slope can be seen by comparing pale blue and turquoise or beige and brown ranges in Fig. 12. At low redshifts, variations with steeper MZR slope lead to ∼10 times higher SFRD at low metallicity. This is due to the presence of a more extended low metallicity tail of the f_SFR(Z,z) (see Section 5.1.1). Lower a_MZR lead to more compact f_SFR(Z,z) – in those cases, for the high metallicity variations even the low-mass galaxies contribute to SFRD at Z_O/H above the considered threshold of 10 per cent solar metallicity. This difference is reduced at higher redshifts due to steeper metallicity evolution of variations with lower a_MZR values – the slope of the low metallicity SFH redshift dependence at z ≲ 2 depends on the z ∼ 0a_MZR. The slope of the dependence is also affected by the contribution of starbursts (compare pale blue and beige or turquoise and brown ranges). As expected, if the fraction of starbursts increases with redshift (B18/C17 implementation), the low metallicity cut of the SFH increases faster towards high redshifts. At z ≲0.2 (within the last ∼2 Gyr) there is little difference in the amount of SFRD occurring at low metallicity between the variations with no/low fraction of starbursts and the B18/C17 implementation. The contribution of starbursts only becomes evident at higher redshifts. At the peak of the cosmic SFH variations with the B18/C17 starbursts implementation lead to ≳3 times higher SFRD at Z_O/H < Z_O/H⊙⁻¹ than the corresponding variations with no/low fraction of starbursts As discussed above, at z ≳ 2 the slope of the redshift dependence of the low metallicity SFRD and the location of the peak of this dependence is dominated by the assumptions about α_GSMF. Observational tracers of star formation happening at low metallicity available at high redshifts (e.g. the volumetric rate of long gamma-ray bursts or, potentially, that of black hole – black hole mergers as a function of redshift) could help to discriminate between the different variations of our models. In particular, it could provide information about the contribution of low-mass galaxies to the total SFRD at high redshifts (thereby constraining jointly the low mass SFMR and – to constrain with regular galaxy surveys.

Considering all model variations discussed in this section, the overall spread between the estimates of the SFRD happening below Z_O/H < Z_O/H⊙⁻¹ ranges between ∼1.6 dex around the peak of the cosmic SFH and ∼3 dex at z > 8. At z ∼ 0 there is ∼2.3 dex variation between the different estimates. The uncertainty is dominated by the assumptions about a_MZR. At the peak of the cosmic SFH the spread is somewhat reduced due to steeper metallicity evolution in variations with lower a_MZR. At z ∼ 8 the spread increases to ∼2.5 dex if the low-mass end of the GSMF is allowed to steepen with redshift, and remains close to ∼1.6 dex if only variations with fixed α_GSMF are considered. Note that the overall spread between the low metallicity SFH estimates discussed in this study is significantly larger than ∼0.5–1.2 dex found in Chruslinska & Nelemans (2019) for the same low metallicity threshold. In their study, effectively only the variations due to the uncertain MZR normalization and the redshift evolution of the low mass end of the GSMF (which are also included in our analysis) were explored. This difference highlights the importance of the additional factors discussed in this study for the low metallicity cut of the cosmic SFH. Note that the spread in the low metallicity SFH estimates as discussed above is not only due to differences in the metallicity distribution of the SFRD at different redshifts found for the considered model variations, but also due to their different total SFRD. The latter differences are effectively eliminated in Fig. B1 in Appendix  B, where we additionally show the fraction of the total SFRD that happens at Z_O/H < Z_O/H⊙⁻¹ as a function of redshift for the discussed model variations.

Finally, we note that the low metallicity part of the SFH is very sensitive to the adopted log₁₀(M_*) limit down to which all the empirical relations are extrapolated. This is especially important at z ≲ 1, where the low metallicity star formation is limited to low mass galaxies, and at z ≳ 4, due to the potential steepening of the low-mass end of the GSMF. With the current data there is no guarantee that those relations can be safely extrapolated to galaxy masses below log₁₀(M_*/M_⊙) ≲8. However, given their abundance, completely neglecting the contribution of those galaxies is not a satisfactory solution either. In variations favoring high metallicity (i.e. with high MZR normalization, shallow MZR slope, no/low contribution of starburst galaxies) ignoring the contribution of galaxies with log₁₀(M_*/M_⊙) < 8 would lead to negligible amount of SFRD happening below 10 per cent solar metallicity even beyond the peak of the cosmic SFH. If the SFR–metallicity correlation significantly weakens/disappears at low specific SFR (see additional model variation discussed in Appendix  C and Fig. C4), this may also significantly reduce the estimated amount of SFRD happening below 10 per cent solar metallicity in variations with relatively flat MZR. In those cases, the low/intermediate redshift star formation can only take place below 10 per cent solar metallicity in low mass galaxies that are MZR outliers (see e.g. left bottom panel in Fig. 4 and in Fig. C1), and so the result strongly depends on the adopted FMR extrapolation. Until better constraints on the properties of the low-mass end the star-forming galaxy population are available, those extrapolations are important source of uncertainty for both the total cosmic SFH, as well as for its low metallicity cut.

6.3 A high metallicity cut of the cosmic SFH

The high metallicity cut of the cosmic SFH – i.e. SFRD that occurs at Z_O/H > Z_O/H⊙ as a function of redshift/time for different variations of our model is shown in Fig. 13. This quantity is relatively well constrained in our models when compared to the low metallicity cut of the SFH discussed in the previous section: the overall spread between the estimates obtained with different model variations from this study ranges from ∼0.8 dex at z ∼ 0 to ∼1.8 dex at z ≳7. The biggest difference appears between the estimate from Chruslinska & Nelemans (2019) (orange range) and all the variations from this study, which lead to considerably higher SFRD occurring at high metallicity beyond the peak of the cosmic SFH. As discussed in Section 5.1, this is due to the fact that the non-evolving FMR assumption leads to much slower metallicity evolution when compared to that implied by the extrapolated redshift dependent MZR. The general features that can be seen for all ranges shown in Fig. 13 can be summarized as follows:

• The high metallicity cut of the SFH shows a gradual increase towards the peak of the SFH (z ∼ 2) and a sharp decrease at high redshifts. This reflects the evolution of the total SFRD. However, the increase towards the peak is shallower and the later decrease is sharper – both due to the overall shifting of the SFRD towards lower metallicities at earler epochs.

• The width of the ranges is set by the assumptions about the SFMR and z ∼ 0 MZR: upper edges of all ranges correspond to high metallicity variations, lower edges of all ranges correspond to low metallicity variations. Assumptions about the low mass end of the GSMF have negligible effect on the high metallicity cut of the cosmic SFH.

Comparing the ranges obtained for different assumptions about the MZR slope (brown and beige or pale blue and turquoise), one can see that the results start to deviate at z ≳ 1, where variations with steeper MZR lead to higher SFRD at high metallicity. This is due to slower evolution of the metallicity distribution at z ≲ 4 for those variations, as discussed in Section 5.1.1. Assumptions about starbursts have little effect on the Z_O/H > Z_O/H⊙ metallicity cut of the cosmic SFH (compare light blue/turquoise and beige/brown ranges). This is due to the fact that (given the FMR), those galaxies only contribute to star formation at relatively low metallicities. The fraction of the total SFRD that happens at Z_O/H > Z_O/H⊙ as a function of redshift for the discussed model variations is shown in Fig. B2 in Appendix  B.

Figure 13.

Star formation rate density at high metallicity [above solar, assuming Z_O/H⊙ = 8.83 Grevesse & Sauval (1998)] as a function of lookback time/redshift. The plotted ranges show estimates obtained for different variations of assumptions about the starburst galaxies [blue/turquoise: as in Boco et al. (2021), brown/beige: based on Caputi et al. (2017) – C17 and Bisigello et al. (2018) – B18] and low-mass slope of the redshift z ∼ 0 MZR (a_MZR). The ranges span between the low and high metallicity extremes as defined in the text and assume ∇_FMR0 = 0.27. The orange region shows the corresponding range from Chruslinska & Nelemans (2019). Note that starburst galaxies contribute to SFRD at relatively low metallicity and have little impact on the high metallicity cut of the SFRD.

Open in new tabDownload slide

6.4 The impact of ∇_FMR0 on the low/high metallicity cut of the cosmic SFH

As discussed in Section 3.1.5, the range of extreme reasonable assumptions about ∇_FMR0 and the potential dependence of this parameter for instance, on the metallicity derivation technique, cannot be easily identified with the currently available observational results. We therefore discuss the sensitivity of our result to the choice of this factor. The impact of ∇_FMR0 on the high and low metallicity cut of the cosmic SFH for the example model variations is shown in Figs 14 and 15, respectively. In essence, higher (lower) ∇_FMR0 increase (decrease) the amount of SFRD that happens at low metallicity (and vice versa for the SFRD at high metallicity). Lower ∇_FMR0 means that the difference between the metallicities of highly star-forming galaxies and their quiescent counterparts of the same mass is reduced. As a consequence, if the anticorrelation between the SFR and metallicity is much weaker than with our fiducial assumptions, the impact of starbursts on the low metallicity cut of the SFH is substancially reduced. This can be seen in Fig. 15, where the lower and upper edges of the range obtained with B18/C17 starburst implementation and ∇_FMR0 = 0.17 (dashed lines) roughly match the range obtained for the corresponding model variations assuming low fraction of starbursts (turquoise range). This also means that the impact of the choice of ∇_FMR0 on the variations with lower fraction of starbursts and/or starburst sequence with smaller offset from the SFMR than in B18/C17 implementation, is smaller than in the examples shown in Figs 15 and 14. As discussed in Section 5.1.1, ∇_FMR0 also affects the rate of the overall decrease in metallicity – and so the slope of the different metallicity cuts of the cosmic SFH as a function of redshift. However, while many different factors affect the slope of the low metallicity cut of the SFH below the peak of the cosmic SFH (see Section 6.2), at high redshifts assumptions about the low mass end of the GSMF clearly dominate the uncertainty in the slope regardless of ∇_FMR0.

Figure 14.

Example illustrating the effect of ∇_FMR0 on the high metallicity cut of the cosmic SFH. The plotted range shows estimate obtained with B18/C17 starbursts implementation and a_MZR ∼ 0.3, assuming ∇_FMR0 = 0.27 (same as in Fig. 13). Grey solid (black dashed) lines mark the location of the upper and lower edges of that range obtained with ∇_FMR0 = 0.3 (∇_FMR0 = 0.17).

Open in new tabDownload slide

Figure 15.

Example illustrating the effect of ∇_FMR0 on the low metallicity cut of the cosmic SFH. The brown range shows estimate obtained with B18/C17 starbursts implementation and a_MZR ∼ 0.3, assuming ∇_FMR0 = 0.27 (same as in Fig. 12). Grey solid (black dashed) lines mark the location of the upper and lower edges of that range obtained with ∇_FMR0 = 0.3 (∇_FMR0 = 0.17). Faint lines in the corresponding colours show the evolution of the upper edge under the assumption of non-evolving low-mass end of the GSMF. For comparison, we also show the range obtained with Boco et al. (2021) starbursts implementation, a_MZR ∼ 0.3 and ∇_FMR0 = 0.27 (turquoise, same as in Fig. 12). If the SFR–metallicity correlation is much weaker than implied by our fiducial ∇_FMR0 = 0.27, the impact of starbursts on the low metallicity cut of the SFH is substantially reduced.

Open in new tabDownload slide

7 CONCLUSIONS

In this study we expand the framework introduced in Chruslinska & Nelemans (2019) to construct the observation-based distribution of the cosmic SFRD at different metallicities and redshifts. We evaluate its uncertainty in light of additional factors that were not discussed in the original study: the redshift-invariant mass–metallicity-SFR relation (FMR), starburst galaxies and the slope of the local mass–metallicity relation. We show that all of those factors have a strong impact on the low metallicity cut of the cosmic SFH – crucial in the discussion of the origin of various energetic transients. Similarly to Chruslinska & Nelemans (2019), we consider different MZR and SFMR variations to cover the range of possibilities discussed in the literature and explore the realistic extremes of the f_SFR(Z,z). The diversity of observational derivations of the MZR and SFMR stems from differences in the applied metallicity and SFR determination techniques – those differences in the methods also affect the FMR determinations. To ensue the consistency of our method, we develop a phenomenological model to construct the FMR based on the local MZR and SFMR. We introduce a parameter ∇_FMR0 to characterize the strength of the SFR–metallicity correlation in the high SFR part of the FMR. That is the only parameter which is not constrained by the choice of z ∼ 0 MZR and SFMR within our description and we discuss the sensitivity of our results to this factor. Our description captures all characteristic features of the observationally inferred relation. It reproduces the known systematic differences between the FMR derived with various metallicity and SFR determination techniques and predicts dependence of the properties of the FMR on the parameters of the z ∼ 0 MZR and SFMR that can be verified with future observational studies. Furthermore, we use the results of Caputi et al. (2017) and Bisigello et al. (2018) to describe the properties of starburst galaxies and include their contribution in our f_SFR(Z,z) models. We contrast the f_SFR(Z,z) derived that way with the f_SFR(Z,z) obtained with the description of starbursts as recently used in Boco et al. (2021) (based on Sargent et al. (2012) and Béthermin et al. (2012), which predict a much smaller contribution of strabursts to the total cosmic SFRD budget). In addition to full f_SFR(Z,z) distributions obtained for different variations of our model, we discuss the resulting cosmic SFH and its low and high metallicity cuts. Our main conclusions are the following:

• As recently discussed in Boco et al. (2021), a redshift-invariant FMR can lead to drastically different conclusions about the metallicity at which the star formation occurs at z ≳ 3 than the redshift-dependent MZR extrapolated to high redshifts (as used in Chruslinska & Nelemans (2019)). However, while the high metallicity cut of the SFH at z .5 is strongly affected by this assumption (Fig. 13), the low metallicity SFH estimate from Chruslinska & Nelemans (2019) falls within the range covered by the variations explored in this study (Fig. 12). In that case, the much shallower metallicity evolution implied by the non-evolving FMR is counterbalanced by the low metallicity contribution of starbursts (assuming that they follow the same FMR as regular star forming galaxies).

• The low metallicity tail of the f_SFR(Z,z)f_SFR(Z,z) is the most sensitive to variations in the model assumptions. We find a factor of ∼200 (∼40) difference in the amount of SFRD occurring below 10 per cent solar metallicity at z ∼ 0 (z ∼ 2) between the extreme model variations considered in this study.

• At low redshifts, the uncertainty of the low metallicity SFH is dominated by the assumptions about the slope and normalization of the local MZR, while both MZR and starbursts strongly affect this quantity at z > 1. However, the importance of starbursts for the low metallicity cut of the SFH is highly sensitive to the assumed strength of the SFR–metallicity correlation (∇_FMR0).

• Beyond the peak of the cosmic SFH, the shape of the low metallicity cut of the SFH (peak location, slope) is set by the assumptions about the low-mass end of the GSMF. Those assumptions also dominate the redshift dependence of the total SFRD at z ≳ 3. Given this sensitivity, we propose that a tracer of the SFR happening at low metallicity available at high redshifts (e.g. volumetric rate of long gamma-ray bursts, or potentially – stellar black hole mergers) could provide constraints on the contribution of low mass galaxies – difficult to constrain otherwise – to the SFRD at those redshifts.

• The uncertainty in the total SFRD in our models at z ≲ 1 is dominated by the uncertain shape of the high mass end of the SFMR. Assumptions about the contribution of starbursts dominate the uncertainty of the cosmic SFH at 1 < z ≲ 3. Model variations with SFMR with sharp flattening at high masses, non-evolving low mass slope of the GSMF and no/small fraction of starbursts may underpredict the total SFRD at z ∼3. Model variations which account for the contribution of starbursts with B18/C17 based prescription lead to ∼2.5 times higher SFRD at z ≳4 and shallower redshift evolution of the total SFRD at z ≳2, more in line with the recent SFRD determinations at high redshifts.

ACKNOWLEDGEMENTS

We thank Karina Caputi and Jarle Brinchmann for helpful discussions. We thank the anonymous referee for their thoughtful comments that helped to improve this paper. MC and GN acknowledge support from the Netherlands Organisation for Scientific Research (NWO). AL and LB are supported by PRIN MIUR 2017 prot. 20173ML3WW, ‘Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes’, and by the EU H2020-MSCA-ITN-2019 Project 860744 ‘BiD4BESt: Big Data applications for black hole Evolution STudies.’

DATA AVAILABILITY

The results of our calculations (allowing to reproduce all of the results shown in the paper) for all model variations discussed in this study are publicly available under this url: https://ftp.science.ru.nl/astro/mchruslinska.

Footnotes

REFERENCES

APPENDIX A: FMR VARIATIONS – ADDITIONAL FIGURES

In this section, we show additional figures illustrating the sensitivity of Z_O/H(SFR,M_*) to the parameter β (Fig. A1) and the Z_O/H(SFR,M_*) dependence on various choices of the parameters of our observation-based model. Fig. A2 shows the values of β as a function of log₁₀(SFR) obtained as described in Section 3.1.4 for different choices of the z ∼ 0 MZR, ∇_FMR0, and SFMR used in this study. Figs A3 and A4 show the comparison between the Z_O/H(SFR,M_*) obtained for different SFMR and ∇_FMR0 choices, respectively.

Figure A1.

Dependence on β parameter in equation (1), illustrated using the best fit MZR z∼ 0 parameters from Curti et al. (2020) and various choices of β. If β can be constrained to values β ≳1, the Z_O/H uncertainty induced by the choice of this parameter is relatively small and affects a relatively small range of galaxy masses. However, if β ≲1, small variation in this parameter can add a appreciable uncertainty to Z_O/H.

Open in new tabDownload slide

Figure A2.

Parameter β as a function of log₁₀(SFR) (see equation 1), depending on the choice of z ∼ 0 MZR (different columns), ∇_FMR0 (different rows), and the shape of the high-mass end of the z ∼ 0 SFMR (different lines). Single value of β provides a good description if the SFMR does not significantly depart from a single power law. See Section 3.4 for the description of different model variations.

Open in new tabDownload slide

Figure A3.

Same as Fig. 4, but for different choices of the z ∼ 0 SFMR and MZR. Left-hand (right-hand) panels assume SFMR with no (sharp) flattening at the high-mass end. Top (bottom) panels assume S20 (PP04) z ∼ 0 MZR. All panels assume ∇_FMR0 = 0.27.

Open in new tabDownload slide

Figure A4.

Same as Fig. 4, but for different choices of ∇_FMR0. Top (bottom) panels assume S20 (PP04) z ∼ 0 MZR. All panels assume SFMR with a moderate flattening at high masses (a_SFMR = 0.72 at the high-mass end).

Open in new tabDownload slide

APPENDIX B: COSMIC SFH – ADDITIONAL MATERIALS

Table B1 provides the coefficients of the broken power-law fits to the total cosmic SFH for the variations with non-evolving low-mass end of the GSMF. Table B2 (B3) provides the corresponding coefficients obtained for the low metallicity cut of the cosmic SFH, fitted for the variations with redshift-dependent (non evolving) low-mass end of the GSMF. Fig. B1 (B2) shows the fraction of the total SFRD that happens below 10 per cent solar metallicity, i.e Z_O/H < Z_O/H⊙⁻¹ (above solar metallicity, i.e. Z > Z_O/H) as a function of redshift/lookback time, obtained for the different variations of our model discussed in Section 6.2 (6.3). Note that those variations in general lead to different total SFRD at different redshift, which also affects the spread between the estimates of the low (high) metallicity part of the cosmic SFH shown in Fig. 12 (13). Fig. B1 allows to directly compare the differences in the low (high) metallicity tail of the f_SFR(Z,z) distribution.