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

The question how much star formation is occurring at low metallicity throughout the cosmic history appears crucial for the discussion of the origin of various energetic transients, and possibly - double black hole mergers. We revisit the observation-based distribution of birth metallicities of stars (f$_{\rm SFR}$(Z,z)), focusing on several factors that strongly affect its low metallicity part: (i) the method used to describe the metallicity distribution of galaxies (redshift-dependent mass metallicity relation - MZR, or redshift-invariant fundamental metallicity relation - FMR), (ii) the contribution of starburst galaxies and (iii) the slope of the MZR. We empirically construct the FMR based on the low-redshift scaling relations, which allows us to capture the systematic differences in the relation caused by the choice of metallicity and star formation rate (SFR) determination techniques and discuss the related f$_{\rm SFR}$(Z,z) uncertainty. We indicate factors that dominate the f$_{\rm SFR}$(Z,z) uncertainty in different metallicity and redshift regimes. The low metallicity part of the distribution is poorly constrained even at low redshifts (even a factor of $\sim$200 difference between the model variations) The non-evolving FMR implies a much shallower metallicity evolution than the extrapolated MZR, however, its effect on the low metallicity part of the f$_{\rm SFR}$(Z,z) is counterbalanced by the contribution of starbursts (assuming that they follow the FMR). A non-negligible fraction of starbursts in our model may be necessary to satisfy the recent high-redshift SFR density constraints.


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. Klencki et al. 2018;Giacobbo et al. 2018). This makes the modelled properties of the population of such systems particularly sensitive ★ E-mail: m.chruslinska@astro.ru.nl to the assumed distribution of the cosmic star formation rate density (SFRD) at different metallicities and redshifts, f SFR (Z,z) (e.g. Neĳssel 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 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  (briefly introduced in Sec. 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 Ω =0.3, Ω Λ =0.7 and H 0 =70 km s −1 Mpc −1 , assume a (universal) Kroupa (2001) initial mass function (IMF) and solar metallicity of 12 + log 10 (O/H) = / =8.83 and Z = 0.017 (Grevesse & Sauval 1998).

AN OBSERVATION BASED ( , ) DETERMINATION
We construct f (Z,z) building on the framework detailed in . 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 * . All distributions are redshiftdependent. 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 (Z,z) is obtained by assigning a metallicity to each * . The SFR and metallicity distributions are modelled as log-normal distributions centered 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 ( / =12+log 10 (O/H) 1 ), and that is the metallicity measure used in our study. Additional scatter in / 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 ( ) and * . 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 Sec. 2 in Chruslinska & Nelemans 2019 and references therein) In this study we explore an alternative method to obtain the redshiftdependent 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 Sec. 3 and Sec. 4 respectively.

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 Mannucci et al. 2010). The relation implies the mass-metallicity correlation, as observed in the MZR, but it also introduces an anti-correlation between metallicty and SFR, such that galaxies of the same stellar mass showing higher than average SFR also have lower metallicities (i.e. the 1 In contrast to the total metal mass fraction or iron abundance commonly used in the simulations/theoretical studies. / is often used as a proxy for the total , but note that this requires assuming a particular (typically solar-like) abundance pattern. galaxy's offset from the average mass-metallicity relation -MZR-is anti-correlated with its offset from the average SFR-mass relation). The existence of such SFR-metallicity anti-correlation has been reported in numerous observational studies (e.g. Mannucci et al. 2011;Yates et al. 2012;Andrews & Martini 2013;Bothwell et al. 2013;Lara-López et al. 2010;Salim et al. 2014;Zahid et al. 2014b;Yabe et al. 2015;Hunt et al. 2016;Sanders et al. 2018;Cresci et al. 2019;Curti et al. 2020;Sanders et al. 2020) and established up to redshift ∼ 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é et al. 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 anti-correlation between the galaxy offsets from the MZR and SFMR is expected to hold (e.g. Torrey et al. 2018). Those timescales 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 andSanders et al. 2020 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 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 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. 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 Sec. 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. (2020) show that when the differences in the properties of the ionized gas within HII 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. (2020) is consistent with the non-or weakly evolving FMR within the redshift and mass range probed in their study. The existence of a 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 normalisation).
In light of the above discussion, we conclude that the assumptions about the metallicity evolution made in , 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 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 a 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 is then defined by the local FMR and the SFMR and GSMF (whose 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 ) 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 Sec. 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  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).

Constructing the FMR from ∼ 0 scaling relations
To construct our model FMR, we describe its 2D projection on the Z / -log 10 (M * ) plane. This projection is the most commonly shown in observational studies. We refer to it as / ( , * ) (see Fig.  1 for the illustration). At fixed SFR, the shape of Z / (SFR,M * ) is found to show the same characteristic features as the MZR: the low mass part of the relation is almost linear (Z / ∝ log 10 (M * ), it bends around a certain turnover mass and flattens at high masses, approaching a constant metallicity value (see e.g. Mannucci et al. 2010Mannucci et al. , 2011Yates et al. 2012;Andrews & Martini 2013;Hunt et al. 2016;Telford et al. 2016;Cresci et al. 2019;Curti et al. 2020;Sanders et al. 2020). To describe this dependence, we use the parametrisation from Curti et al. (2020): where / ;0 is the asymptotic metallicity of the high-mass end of the relation, 0; is the turnover mass, is the slope of / ( , * ) at * << 0;SFR and regulates the width of the knee of the relation. The / ( , * ) dependence on SFR is expected primarily in the turnover mass 0;SFR , as illustrated in the right panel of Figure 1. The left panel of Figure  1 illustrates the average ∼ 0 MZR and SFMR. We relate the / ( , * ) parameters to those of the average ∼ 0 scaling relations -each of the parameters of eq. 1 is discussed in the relevant subsection below. Discussion of the observational properties of the / ( , * ) in the reminder of this section and our choices are motivated by the results shown in the references given at the beginning of this section -we do not list them again in each of the subsections, unless we refer to a result from a specific paper(s).
Equation 1 describes the shape of the relation, but the main observable property that characterizes the FMR is the quantity that relates the galaxy's offset from the average MZR and SFMR (i.e. quantifies the strength of the SFR-metallicity correlation). To describe the strength of the SFR-metallicity correlation we introduce the coefficient ∇ FMR 2 , defined as follows: where = Z / , . -Z / , and = log 10 (SFR . ) -log 10 (SFR ) are the galaxy's offsets from the average MZR and SFMR . Higher values of ∇ FMR imply a stronger correlation.

The asymptotic metallicity
A clear feature present among the observational / ( , * ) determinations is the high mass flattening. We note that some authors report a reversal of the (anti-)correlation 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 / ( , * ) and include it within our model. Observationally, the value of the asymptotic metallicity / ;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 ( / ; 0 ). As such, it is affected by the choice of the metallicity determination method, which leads to systematic offsets in the normalisation of the MZR (compare examples shown in Telford et al. 2016 and in

The slope at low masses and high SFR
At the low mass end, the Z / ( , * ) for fixed log 10 (SFR) seems to be well approximated with a linear dependence. Accordingly, in this regime ( * << 0; ), eq. 1 reduces to: The slope of this dependence is generally steeper than that of the average ∼ 0 MZR. seems to be affected by the choice of the metallicity determination technique (see e.g. examples shown in Yates et al. 2012;Hunt et al. 2016;Telford et al. 2016;Cresci et al. 2019), but it is difficult to read off any systematic trend based on the examples shown in the literature (especially keeping in mind the presence of other differences in the methods that also affect the relation). Z / ( , * ) obtained for different log 10 (SFR), show no clear SFR dependence and appear nearly parallel to each other. For simplicity, here we assume that is independent of SFR. However, we note that in their recent study, Curti et al. (2020) find that decreases with decreasing SFR to values similar to that of their average z∼0 MZR slope (see Fig. 9 therein). This can lead to significantly different metallicities of low M * , low SFR galaxies when the FMR is extrapolated down to low masses ( 10 8 M ). In appendix C we introduce a additional variation of our model that accounts for such a dependence and discuss its influence on our results. We can recover the characteristics summarised above by applying eq. 2 in the low-mass regime and assuming that ∇ FMR =∇ FMR0 is fixed in this regime. The slope then follows directly from the slopes of the MZR and SFMR and the correlation coefficient ∇ FMR0 . For Z / = a MZR log 10 (M * ) + b MZR and log 10 (SFR) = a SFMR log 10 (M * ) + b SFMR (appropriate at * << 0;SFR ), we thus obtain: The slope of the / ( , * ) relation at fixed SFR is then related to the slopes of the low mass parts of the MZR ( MZR ) and SFMR ( ) and the parameter describing the strength of the SFR-metallicity correlation ∇ FMR0 : The values MZR and are set by the choice of the local scaling relations. The choice of ∇ FMR0 is discussed in Sec. 3.1.5. The dependence of on those parameters in the relevant range of values is shown in the right top panel of Fig. 2.
can be expressed with the parameters of the z∼0 average MZR as:

Turnover mass as a function of SFR
Observational studies suggest that 0;SFR increases with SFR (i.e. for higher SFR, / ( , * ) flattens at higher masses). The steepness of this dependence governs the spacing between the Z / (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 0; ). Curti et al. (2020) find that a linear dependence (log 10 ( 0; )∝ log 10 (SFR)) can well describe the trend seen for their z∼0 galaxy sample. We note that the linear dependence between log( 0;SFR ) and log 10 (SFR) is a natural consequence of the assumptions listed earlier in this section, that and / ;0 are independent of SFR. In the * << 0; regime eq. 1 reduces to eq. 3, which in our The dependence of on MZR , SFMR and ∇ FMR0 is shown in Fig.  2.

The value of
describes the bending of the / ( , * ) relation. The higher the value of , the sharper the transition from the linear to flat regime and the smaller the range of log 10 (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 / ( , * ) relation for different SFR bins. While their best fit relations give different values in different SFR bins, the authors note that there is no clear dependence 4 Within our framework the value of the slope of this dependence is also the value of the parameter that minimizes the scatter in the 2D projection of the FMR that is conventionally used to discuss the strength of the three parameter dependence: / − where = ( * ) − ( ), as originally proposed by Mannucci et al. (2010). Although, as can be seen in eq. 6, is not a good measure of the strength of the correlation between the SFR and metallicity, as it depends strongly on the shape of the two relations, in particular on the slope of the low mass end of the MZR. What really describes the strength of the correlation is ∇ FMR . 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 / ( , * ). This condition is only applicable in the SFR range that is covered by the z∼0 SFMR ( 10 −3 M /yr and 10 2 M /yr). Numerous solutions are allowed -especially at low SFR 10 −1 M / , where the intersection with the z∼0 MZR is reached at log( * )<log( 0;SFR ) i.e. where / ( , * ) is virtually insensitive to the choice of . We use the numerical solutions as a guide and adopt a simplified dependence on log 10 (SFR), verifying that / ( , * ) for galaxies at z∼0 SFMR leads to / 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 / ( , * ) 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.

The strength of the SFR-metallicity correlation
Observationally, ∇ FMR appears to be a function of * 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
/ ;0 appears to be roughly independent of the SFR means that the observed / ( , * ) 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 / correlation in the high * regime (i.e. small/zero ∇ FMR ). Note that, by construction (fixing the value of / ;0 ), we recover such high * weakening of the correlation in our method. A similar behaviour is also seen in the Z / -log 10 (SFR) projection of the FMR, where at fixed log 10 (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 −10.5 /yr). When projected onto the Z / -log 10 (M * ) plane, this means that the correlation is weaker 'above' than 'below' the average ∼0 MZR. That is also the case within our description: for the same absolute value of , ∇ FMR is smaller above than below the ∼ 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 * and low SFR/low sSFR weaker correlation regimes, there is no agreement on the precise region of the FMR (in terms of the Z / , log 10 (M * ) and log 10 (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 10 (M * ) and sSFR depending on the parameters of the ∼ 0 MZR and SFMR. We stress that we only assume that ∇ FMR is fixed (i.e. ∇ FMR =∇ FMR0 , independent of * and SFR) in the low mass and high SFR part of the FMR -as guided by the observed characteristics of / ( , * ) (but see Appendix C, where we further relax this assumption and introduce a 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 / ( , * ) 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), Salim et al. (2014) find ∇ ∼0.3 in the strong correlation regime. Applying a infrared based SFR derivation method, they find somewhat lower values of ∇ ∼0.2. In their study focusing on ∼ 2.3, Sanders et al. (2018) use several example metallicity calibrations (both strong line-based/theoretical and 'direct'/empirical) and find ∇ ∼0.11 -0.27, hinting at potential dependence on the choice of metallicity indicator and calibration. However, as argued by Sanders et al. (2020), 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. (2020), where different empirical metallicity calibrations are used at z∼0 and at higher redshifts. They find best-fit ∇ FMR =0.27 at z∼0 and somewhat shallower dependence at z∼2.3 (∇ ∼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 10 (M * ) and redshift bins as found in the IllustrisTNG simulations (see Table 1 therein). They report values in the range ∇ ∼0.25 -0.34 (except for the lowest log 10 (M * )∼9 and highest mass bin log 10 (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. (2020) as our fiducial choice. Given the uncertainty of this parameter, we consider values between 0.17 -0.3 dex to discuss the sensitivity of our results to this choice 5 .
3.1.6 Calculating metallicities of galaxies at different with / ( , * ) / ( , * ) constructed as described in this Section is fully determined by the choice of ∇ FMR0 and the parameters of the ∼ 0 MZR and SFMR. We introduce the variations of those local relations that we explore in this study in Sec. 3.4. We further assume that / ( , * ) 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 centered around the redshift-dependent SFMR, as discussed in Sec. 2. We then use / ( , * ) to assign the metallicity. Furthermore, observational studies indicate the residual scatter around the FMR ∼0.05 dex. To account for that, we add a normally distributed scatter =0.05 dex to metallicities assigned with our / ( , * ).

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 anti-correlated, i.e. setting the coefficient from eq. 2 to ∇ = MZR SFMR ≈ 0.33, where MZR =0.1 dex and SFMR =0.3 dex describe the scatter around the average MZR and SFMR respectively. The same ∇ FMR is used independent of * , SFR or redshift. In the strong correlation regime discussed in Sec. 3.1.5 (at low/intermediate * and high SFR) and at ∼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 * and low SFR/sSFR. This behaviour is not captured with the simple description given by eq. 2 with fixed ∇ FMR . Another important difference is that, in this study the Z / (M * ,SFR) (once defined with ∼ 0 relations) is assumed to be redshift-invariant. Chruslinska & Nelemans (2019) use eq. 2 to calculate offsets relative to SFMR and MZR as found at any given redshift. This means that, due to strong MZR evolution at 2, their Z / (M * ,SFR) is redshift-dependent.

FMR example: comparison with Sanders et al. (2020)
In this Section we present an example / ( , * ) obtained with our phenomenological model, where we keep all the input assumptions as similar as possible to Sanders et al. (2020). We choose 5 Lower values of ∇ FMR were also reported, but we find that ∇ FMR0 0.17 underpredict the redshift evolution of the MZR at 1.5 (where the impact of biases discussed in Sec. 3 is not so severe, and so its evolution is reasonably constrained). Much lower ∇ FMR0 also make it difficult to reproduce the width of the ∼ 0 MZR of ∼0.1 dex (unless the scatter around the FMR is in reality higher than the typically reported ∼0.05 dex). this study, as the authors provide estimates of all the ingredients that are necessary to construct our / ( , * ) (the ∼ 0 MZR, SFMR and ∇ FMR0 ) 6 , which allows for a direct comparison with their results. The resulting / ( , * ) is shown in Fig. 3. In the bottom panels we include two other common 2D projection of the FMR: the / -log 10 (SFR) plane (bottom left panel) and the / -log 10 (sSFR) plane (bottom middle panel). The thick gray line in the main panel of Fig. 3 shows the ∼ 0 MZR as given in Sanders et al. (2020). The coloured points around that line represent the ∼ 0 population of star forming galaxies from our model 7 . The black line indicates a fit to the maximum density region occupied by the ∼ 0 galaxy sample described above, to ensue that the ∼ 0 MZR is reproduced. The right bottom panel shows the / residuals around the MZR obtained that way for several example log 10 (M * ), which shows that the typically indicated intrinsic width of the relation ∼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 / ( , * ) for various fixed log 10 (SFR) values. The faint coloured solid lines in the background were obtained with the best-fit FMR given by Sanders et al. (2020) (see eq. 10 therein), plotted roughly in the range of * and SFR probed by their galaxy 7 * are sampled from our ∼ 0 GSMF, SFR are sampled from SFMR from Sanders et al. (2020) with Gaussian scatter and the metallicity is assigned using our / ( , * ) with scatter Figure 3. Example / ( , * ) obtained with our method (coloured lines -different log 10 (SFR) values) and using the ∼ 0 MZR, SFMR and ∇ FMR0 from Sanders et al. (2020) (S20). Colored 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 gray line and shaded area -3 MZR region) is well reproduced with our / ( , * ) and z∼0 GSMF and SFMR. The gray dashed lines show fits to the maximum density of galaxies sampled from GSMF and SFMR at =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 / ( , * ) as would be obtained with eq. 2 and ∇ FMR0 =0.27 used at all masses and SFRs. Projections onto / -log 10 (SFR) and / -log 10 (sSFR) planes are shown in the bottom panels for several log 10 (M * ) values (black lines -log 10 (M * ) from 6 to 12). In the bottom left panel we indicate the SFR corresponding to SFMR value for a given M * at = 0 (cross), 2.3 (square) and 3.3 (triangle). Red dots indicate the SFR of a galaxy ±3 away from the SFMR at ∼ 0. The rightmost bottom panel shows the offsets of ∼ 0 model galaxies from the ∼ 0 MZR, plotted for several mass bins. sample. We also plot their ∼ 2.3 (squares) and ∼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 10 (SFR) within the uncertainties provided by the authors for each of the stacks. The long dashed and dotted gray lines indicate the projected ∼ 2.3 and ∼ 3.3 MZR, calculated in the same way as the black line, but using the ∼ 2.3 and ∼ 3.3 SFMR from Sanders et al. (2020) respectively, and our GSMF estimated at the corresponding redshifts. The overall agreement is remarkable. We also show / ( , * ) for various fixed log 10 (SFR) values as would have been obtained with the approach used in  i.e., using eq. 2 with ∇ FMR =∇ FMR0 =0.27 fixed (i.e. the same at all SFRs and masses, see colored dashed lines in the main panel of Fig. 3). It can be seen that the / ( , * ) constructed that way is identical with our model at low/intermediate * and high SFRs (the strong correlation regime), and starts to deviate at high * and low SFR (above the ∼0 MZR), where ∇ < ∇ 0 . This demonstrates that the mass and SFR dependence of the SFR-metallicity correlation discussed in Sec. 3.1.5 is present in our description.

Model variations
In this study we consider several variations of the base ∼ 0 MZR and SFMR relations, representing extreme choices of the shapes of the two relations and the MZR normalisation  , * ) at fixed log 10 (SFR) (increasing from left to right, see legend). Different panels correspond to different choices of the ∼ 0 MZR (indicated by the thick gray line). All panels assume SFMR with a moderate flattening ( ∼ 0 SFMR with a =0.72 at high masses and a =0.83 at low masses) and ∇ FMR0 =0.27. Gray dashed lines indicate the projected MZR evolution as obtained by performing fits to maximum density of galaxies in 12+log 10 (O/H) -log 10 (M * ) plane, where the galaxies were sampled from GSMF (assuming = ) and SFMR at several higher redshifts and their Z / was assigned using the / ( , * ) shown in each of the panels. Horizontal orange lines indicate solar and 10% solar metallicity, assuming Z / =8.83 (Grevesse & Sauval 1998 2019 for the details). This allows us to explore the extremes of the f (Z,z) distribution. We briefly introduce each of the MZR and SFMR variations used to construct our Z / (SFR,M * ) below.

MZR variations
The MZR is parametrised in the same way as Z / (M * ,SFR), i.e.: MZR parameters for the variations considered in this study are given in the top rows of between the metallicities derived with the theoretical calibrations (e.g. Kobulnicky & Kewley 2004) and the so-called direct method or empirical calibrations (e.g. Pettini & Pagel 2004;Sanders et al. 2020), where the latter yield estimates that are ∼2 times lower than the former. The choice of KK04 and PP04 calibrations maximises the difference in / ; 0 (and in the low/high metallicity tails of the f (Z,z), see Chruslinska & Nelemans 2019), i.e. metallicity estimates obtained with other methods fall in between (see e.g. Fig  15 in Maiolino & Mannucci 2019) 8 . Furthermore, we consider two additional variations: S20 and KK04b, that have / ; 0 as in the PP04 and KK04a variations respectively but a shallower low mass slope ∼ 0.3. The S20 variation is based on the recent direct-method based determination from Sanders et al. (2020). The parameters of the KK04b variation were chosen in such a way that the MZR is identical with the one in KK04a variation at log 10 (M * / M ) 8.7 and differs at lower masses (where the original determination was not constrained by the data). Our motivation for including those variations is the following: the shape of the z∼0 MZR, in particular its low mass end slope and log 10 0; 0 is what sets and log 10 M 0;SFR within our framework, as the dependence on the SFMR parameters is reduced by multiplication by ∇ 0.3. However, the MZR slope and turnover mass are not well constrained with the current data (the purely linear regime is often not probed) and depend on the adopted parametrisation, for the same galaxy sample one can fit different a and turnover masses (e.g. Curti et al. 2020, for the parametrisation used in eq. 1, is also somewhat degenerate with the MZR turnover mass and slope). Recent MZR determinations by Curti et al. (2020) and Sanders et al. (2020) both indicate shallower a ∼0.3 than typically found in earlier studies (a ∼0.6, e.g. Andrews & Martini 2013;Zahid et al. 2014b). Curti et al. (2020) and Sanders et al. (2020) discuss several potential biases in earlier studies that may be causing this difference (e.g. relatively bright, high SFR galaxies may dominate the low mass end of the sample and therefore induce a bias towards lower metallicity). Alternatively, massive star-based metallicity determination methods (yielding / ; 0 consistent with empirical/direct methods) lead to a ∼0.6, similar to earlier studies (see discussion in Sanders et al. 2020). Given those differing results and until better observational constraints on the low mass end of the MZR are available, we consider 0.3 a 0.6 as realistic. Note that this uncertainty in the MZR slope at M * 10 8.5 M leads to significant differences in metallicity when the relation is extrapolated down to low stellar masses (e.g. at log 10 (M * / M )=6 there is ≈0.56 dex difference between the metallicity estimated with KK04a and KK04b z∼0 MZR variations and ≈0.76 dex difference if S20 and PP04 z∼0 MZR variations are compared). By considering a range of slopes we explore the uncertainty associated with the low mass MZR extrapolation on the (low metallicity part of) f SFR (Z,z) estimate. In this study, the MZR is only used to construct / ( , * ) and we do not need to describe its evolution with redshift.

SFMR variations
To describe the ∼ 0 SFMR we follow  and assume the low/intermediate mass slope and normalisation based on Boogaard et al. (2018). As discussed in , 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 : 'no flattening' -a linear relation between log 10 (SFR) and log 10 (M * ) with a single slope at all masses and 'sharp flattening'with a ∼ 0 at high masses (using the parametrisation from Tomczak et al. (2016) 9 ). 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 to 2 < at log 10 (M * )=log 10 M 0;SFMR0 . Variation moderate/S20 assumes a high mass slope as in the ∼ 0 SFMR shown in Sanders et al. (2020), 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 .

The FMR for the considered model variations
In this section we summarize the key differences between / ( , * ) obtained for the different variations of the ∼0 MZR, SFMR and ∇ FMR0 considered in this study. It can be seen from eq. 6 that the impact of the SFMR parameters on the SFR-dependent turnover mass log 10 (M 0;SFR ) (and on ) is reduced by the multiplication by ∇ 0.3. Considering a range of ∼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 =0.83 would only affect by 0.08. This is illustrated by the gray bands in the top panels in Fig. 2. At the same time, for the range of 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 / ( , * ). ∇ FMR0 also visibly affects the relation, while SFMR has a relatively mild impact on its shape within our framework. / ( , * ) obtained for the different variations of the ∼ 0 MZR considered in this study and for our fiducial ∇ FMR0 =0.27 is shown in Fig. 4. The shift in normalisation 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 10 (SFR) lines depending on MZR (compare right and left 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 / ( , * ) variation with a 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 10 (SFR) relation for different choices of the local scaling relations). In practice, cases with a 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 9 Parametrised as follows: log 10 (SFR) = s 0 − log 10 1 + ( M * M 0;SFMR0 ) −a SFMR weakly affected by the choice of ∇ FMR0 (shifting towards smaller values with decreasing ∇ FMR0 ). We note that using the ∼ 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.

STARBURST GALAXIES
A common approach to describe the SFR distribution of galaxies is to use a gaussian distribution centered around the redshift-dependent SFMR 10 . 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 11 than those of the regular star forming galaxies of the same mass and redshift. Several authors estimate the fraction of starburst galaxies (f ; 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 ∼2-3% (e.g. Rodighiero et al. 2011;Sargent et al. 2012;Béthermin et al. 2012;Ilbert et al. 2015;Schreiber et al. 2015). Despite the relatively low f values, the starburst's contribution to the total cosmic SFRD could still amount to ∼ 10% 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  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 determinations are based on galaxy samples limited to relatively massive objects (log 10 (M * / M )>10) and 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 10 (SFR) -log 10 (M * ) plane to much lower masses and higher redshifts. While their results are consistent with the previous determinations at high log 10 (M * ), they show that f is a strong function of stellar mass (increasing towards lower log 10 (M * )) and (increasing with redshift). Moreover, they indicate that starburst galaxies follow a distinct sequence in the log 10 (SFR)log 10 (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. Sargent et al. 2012;Béthermin 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. 10 But see Boco et al. 2019Boco et al. , 2021 for alternative method based on the galaxy SFR functions (number density of galaxies per logarithmic bin of SFR). 11 Note that various criteria are used in the literature to distinguish starbursts and regular star forming galaxies. Most commonly the criteria are based on the SFR and require that the starburst SFR is at least a factor of a few (factors between 3 and 10 are used in the literature) higher than that of an average galaxy of the same mass and at the same redshift (e.g. Orlitova 2020) Therefore, they affect the low metallicity tail of the f (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 (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.

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 10 (M * ) of star forming galaxies from the galaxy stellar mass function as described in .
• We use f to describe the fractions of starburst galaxies and regular star forming galaxies at each log 10 (M * ) and . Our choice is outlined further in this section.
• The SFR of regular galaxies is given by the SFMR from . The SFR of starburst galaxies at each log 10 (M * ) follows a normal distribution with scatter . The peak of the distribution can be related to the galaxy's log 10 (M * ) with a linear relation: log 10 (SFR SB ) = a SB log 10 (M * ) + b SB , to which we refer as the starburst sequence We set the value of by defining the offset of the starburst sequence from the SFMR, i.e.
= + Δ − . 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 Sec. 3.1).
We consider two sets of parameters describing the properties of starbursts. Firstly, 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 =0.03 (independent of mass and redshift), a starburst sequence with scatter =0.24, located Δ − =0.59 dex above the SFMR and parallel to the SFMR (specifically, we assume the starburst sequence slope of a =0.83, corresponding to the low mass slope of the SFMR from Chruslinska & Nelemans 2019 12 ) In light of the results of Caputi et al. (2017) and Bisigello et al. (2018) this implementation severely underestimates both f 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 and the properties of the starburst sequence. We provide the details of this implementation in Sec. 4.1.1 and 4.1.2 and refer to it as the B18/C17 implementation in the reminder of this paper.

The fraction of starbursts
To describe the mass and redshift dependence of , we use the observational estimates obtained by Bisigello et al. (2018) (provied for three redshift bins: 0.5 < < 1, 1 < < 2, 2 < < 3; see   Table. 2) for different redshift bins. Open symbols mark estimates where the galaxy sample from Bisigello et al. (2018) is below 90% mass completeness.
The assumed fraction of starbursts versus stellar mass and redshift is shown in Fig. 5. At each redshift bin covered by the data, we assume the following relation between and log 10 (M * ): All stellar masses in eq. 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 at 0.5 < < 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% at all masses at z=0. At z=0.75,1.5 and 2.5, log 10 (M fix ) is the lower edge of the lowest mass bin where the galaxy sample from Bisigello et al. (2018) is within 90% stellar-mass completeness. We make a conservative assumption and use a fixed value below that mass (i.e., = a 1 log 10 (M fix ) + b 1 ). At = 4.4, log 10 (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, is fixed at the same mass as at z=2.5 (i.e. log 10 (M * )=8.24) and we assume a linear relation to describe the dependence between log 10 (M * )=8.24 and log 10 (M * )=9.25 . This exception is necessary to avoid a decrease in the fraction of starbursts at log 10 (M * )<9.25 at z>2.5, which would break the trend seen in the data (see Figure 5). At log 10 (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 =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 -log 10 (M * ) dependence is well described with the same slope across all four redshift bins. Therefore, for simplicity we  assume 1 = −0.3. The slope of the relation in the higher mass part ( 2 ) and the mass log 10 (M 0; ) separating the high and low mass parts of the relation increase with redshift.

The starburst sequence
We guide our description of the SFR distribution of starburst galaxies with the results shown in Fig. 3 and Fig. 7 from Caputi et al. (2017) and Fig. 6 and Fig. 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 =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 evolution with redshift. For simplicity, we assume =0.94 (average between the best fit values in 3 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 Δ − we focus on the results for the intermediate masses log 10 (M 10 / 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 Δ − 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.

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 (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 (Z,z) cases in terms of the amount of SFRD occurring at low (below 10% solar metallicity; Z / Z / -1) and high (above solar Z / Z / ) 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 (Z,z) cases -obtained with ∼ 0 MZR with low normalisation (PP04 or S20) and SFMR with sharp flattening at high masses, and the 'high metallicity' f (Z,z) cases -obtained with ∼ 0 MZR with high normalisation (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 (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 * 10 8 M ). The importance of the assumed z 3 extrapolation is discussed in Sec. 5.1, where we compare the f (Z,z) obtained with the non-evolving FMR and with the redshift-dependent MZR based approach used in . Metallicity of * 10 8 M galaxies assigned with our FMR is primarily sensitive to the MZR slope and normalisation and the strength of the SFR-metallicity anti-correlation at low masses/SFRs. Those factors are discussed in Sec. 5.1.1 and further in Appendix C. In Sec. 5.2 we discuss the f (Z,z) under different assumptions about the contribution of starburst galaxies.

f (Z,z) with redshift-invariant FMR
The comparison between the f (Z,z) constructed with the use of redshift-invariant Z / (SFR,M * ) and the corresponding distributions from  is shown in Fig.  7. We compare the variations with ∼ 0 MZR and SFMR choices as in the high (left) and low (right) metallicity extremes from Chruslinska & Nelemans (2019) (see Sec. 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 . Fig. 7 shows that the metallicity distributions start to deviate around 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 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 Sec. 3, the extrapolated MZR evolution as assumed by  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 (Z,z) resulting from the extrapolated evolution of the galaxy metallicity distribution with redshift. The extrapolated MZR evolution leads to a ∼ 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 (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 (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 Sec. 4.1 therein), this assumption has relatively little effect on f (Z,z) at 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 (Z,z) is analogous for other model variations.

The effect of the ∼ 0 MZR slope and ∇ FMR0
The effect of the choice of the ∼ 0 MZR slope ( MZR ) and ∇ FMR0 on f (Z,z) is shown in Fig. 9. Similarly to Fig. 7, the left and right panels show high and low metallicity model variations respectively. Top (bottom) panels show variations with ∼0.3 ( ∼0.6). Different contours illustrate the impact of ∇ FMR0 As can be expected, f (Z,z) for variations with steeper ∼ 0 MZR (higher MZR values) have a more extended low metallicity tail at each redshift. However, variations with lower MZR values feature a much steeper metallicity evolution at 4. At higher redshifts the metallicity peak of f (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 MZR . 13 This can be understood by looking at the comparison of Z / (SFR,M * ) obtained for different ∼ 0 MZR, shown in Fig. 4. In this plane, galaxies with the same log 10 (M * ) at higher redshifts shift to higher log 10 (SFR) lines. The rate of this shifting 13 For instance, for the ∇ 0 = 0.27 cases shown in Fig. 9, the peak metallicity decreases by ∼0.55 dex (∼0.4 dex) between z=0 and z=4 in the model with S20 (KK04b) ∼ 0 MZR. This decrease is only ∼ 0.1 dex (∼0.15dex) for the models with steeper PP04 (KK04a) ∼ 0 MZR. Between z=4 and z=10, the peak metallicity decreases by ∼ 0.1 dex (∼0.18 dex) for the cases with SFMR with no (sharp) flattening (irrespective of the MZR ).  (2019). All models assume GSMF with non-evolving low mass end. Orange horizontal dashed lines indicate solar, 10% solar and 1% solar metallicity (assuming Grevesse & Sauval 1998 solar metallicity scale). The right metallicity axis was obtained from / 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.
is dictated by the redshift-dependent SFMR. The offset in Z / between the different log 10 (SFR) curves at fixed log 10 (M * ) is bigger for variations with lower MZR (compare left and right panels in Fig 4) and similarly for higher ∇ FMR0 (see Fig. A4 in the Appendix A). This translates into steeper decrease of the f (Z,z) metallicity peak with redshift in the variations with lower MZR and higher ∇ FMR0 . At sufficiently high redshifts (where almost all galaxies occupy the high SFR, linear regime of Z / (SFR,M * )) the spacing in Z / for different log 10 (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 4 and fixed log 10 (M * ) is essentially the same. However, the full f (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 a SFR-dependent ∇ FMR in which the SFR-metallicity correlation is assumed to disappear (i.e. the FMR breaks down) at SFR corresponding to ∼ 0 SFMR galaxies with * 10 8 M . Overall, this assumption has minor impact on the estimated f SFR (Z,z) distribution (see Fig. C2).

f (Z,z): the impact of starbursts
The comparison of f (Z,z) obtained with and without including starbursts is shown Fig. 10. Top and middle panels show the same ∼ 0 MZR and SFMR variations as in Fig. 9. Different contours correspond to different starburst implementations (introduced in Sec. 4.1). The effect of including starbursts in our calculations is twofold: first, a fraction 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). Second, assuming that starbursts follow the FMR, they contribute to star formation happening at relatively low (for their mass and redshift) metallicity. Therefore, including starbursts shifts the peak of the f (Z,z) to lower metallicities and broadens the low metallicity part of the distribution.
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 . The total SFRD is about 2.5 times higher at 4 in the case with starbursts (brown line, bottom panel in Fig. 10). Note that we fix f beyond = 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 . f SFR (Z,z) significantly deviate at 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.
difference with respect to cases without starbursts is negligible. This is expected, given the low fixed f =3% and the fact that in this prescription the starburst sequence is within the 2 scatter of the SFMR. The difference is slightly more pronounced if SFMR with a sharp flattening at high masses is used (right 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).

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 (Z,z) distributions shown in Sec. 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 (Z,z) distribution may lead to considerable uncertainty in the SFH cut for a particular choice of metallicity threshold. The 10% solar/solar Z / 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 drives the uncertainty of the f (Z,z) in different regimes (i.e. low/high metallicity, low/high redshift). As discussed in Sec. 5, the difference between the f (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 is negligible. Therefore, in this section we only discuss the former.

Cosmic SFH
The cosmic SFH obtained for different model variations considered in this study is shown in Figure 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 ( ∼ 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 Figure 11 with a broken power law: SFRD = (1+z) , fitting 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 4, before the sharp decrease at > 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 10 (M * ) < 8 (see discussion in Chruslinska & Nelemans 2019).
• All variations show a sharp decrease in the SFRD at ∼ 7. This is due to the rapid evolution of the GSMF normalisation between ∼ 7 and 8, seen when observational GSMF estimates from different redshifts are combined (see Fig. 3 in Chruslinska & Nelemans 2019).
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 Sec.

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 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 Figure 11 we contrast the total SFRD that results from various realisations 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 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 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. that assume SFMR with no flattening at high masses (and C17/B18 starburst implementation) and the commonly used, best-fit estimate from Madau & Dickinson (2014) (gray 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 < < 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 + ) −2.9 ). This is striking for the upper edges of our estimates (see also Tab. 3). Shallower SFRD evolution than estimated by Madau & Dickinson (2014) is in line with the most recent high redshift observational estimates (e.g. Khusanova et al. 2021;Gruppioni et al. 2020;Loiacono et al. 2021, although the uncertainties in those measurements are still substantial), indicating that a significant fraction of the SFRD at 2 still occurs in dust obscured galaxies (and hence was missed by the earlier -mostly UV light-based surveys). Those estimates 14 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.  (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.
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 10 (M * / M ) = 6 -as discussed above, the result of this extrapolation at high 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 ∼ 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 normalisation 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.

Note on the low metallicity threshold
When referring to low metallicity, throughout this paper we use a limit of <10% solar oxygen abundance (with / 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% solar (depending on the considered model of binary/stellar evolution, e.g. , Broekgaarden et al. (2021. In that sense, depending on the model, the interesting low metallicity threshold could be as high as 50% 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 Table 3. Coefficients of the power law fits to the cosmic SFH (lower and upper edges of the ranges shown in Fig. 11) Sanders et al. 2020). In general, interstellar medium (ISM) of galaxies is enriched in oxygen and iron on different timescales 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 et al. 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 super-solar O/Fe 2 (O/Fe)| in the ISM of ∼ 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 et al. 2004;Tolstoy et al. 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.

Results
SFRD that occurs at Z / <Z / -1 as a function of redshift/time for different variations of our model is shown in Figure 12. Note that the assumptions about the high mass end of the SFMR have no impact on the low metallicity tail of f (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 Figure 12 fitted in several ranges. Those are listed in Tab. B3 and Tab. B2 in the Appendix B. The following features are common to all model variations: • At 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 ∼ 0 MZR: lower/upper edges of all ranges correspond to variations with high/low MZR normalisation (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 normalisation and non-evolving/steepening GSMF low mass end.
• If 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 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.
The impact of the assumed ∼ 0 MZR slope can be seen by comparing pale blue and turquoise or beige and brown ranges in Figure 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 (Z,z) (see Sec. 5.1.1). Lower MZR lead to more compact f (Z,z) -in those cases, for the high metallicity variations even the low mass galaxies contribute to SFRD at Z / above the considered threshold of 10% solar metallicity. This difference is reduced at higher redshifts due to steeper metallicity evolution of variations with lower MZR values -the slope of the low metallicity SFH redshift dependence at 2 depends on the ∼ 0 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 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 / <Z / -1 than the corresponding variations with no/low fraction of starbursts As discussed above, at 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 . 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 GSMF) -difficult 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 / <Z / -1 ranges between ∼1.6 dex around the peak of the cosmic SFH and ∼3 dex at >8. At ∼0 there is ∼2.3 dex variation between the different estimates. The uncertainty is dominated by the assumptions about MZR . At the peak of the cosmic SFH the spread is somewhat reduced due to steeper metallicity evolution in variations with lower MZR . At ∼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  for the same low metallicity threshold. In their study, effectively only the variations due to the uncertain MZR normalisation 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 / < / -1 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 10 (M * ) limit down to which all the empirical relations are extrapolated. This is especially important at 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 10 (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 normalisation, shallow MZR slope, no/low contribution of starburst galaxies) ignoring the contribution of galaxies with log 10 (M * / M )<8 would lead to negligible amount of SFRD happening below 10% 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% solar metallicity in variations with relatively flat MZR. In those cases, the low/intermediate redshift star formation can only take place below 10% 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.

A high metallicity cut of the cosmic SFH
The high metallicity cut of the cosmic SFH -i.e. SFRD that occurs at Z / >Z / 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 ∼ 0 to ∼1.8 dex at 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 Sec. 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 summarised as follows: • The high metallicity cut of the SFH shows a gradual increase towards the peak of the SFH ( ∼ 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 ∼ 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 1, where variations with steeper MZR lead to higher SFRD at high metallicity. This is due to slower evolution of the metallicity distribution at 4 for those variations, as discussed in Sec. 5.1.1. Assumptions about starbursts have little effect on the Z / >Z / 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 / > / as a function of redshift for the discussed model variations is shown in Fig. B2 in Appendix B.
6.4 The impact of ∇ FMR0 on the low/high metallicity cut of the cosmic SFH As discussed in Sec. 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 Fig. 14 and Fig. 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 anti-correlation 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 Fig. 15 and 14. As discussed in Sec. 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 Sec. 6.2), at high redshifts assumptions about the low mass end of the GSMF clearly dominate the uncertainty in the slope regardless of ∇ FMR0 .

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 massmetallicity 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 , 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 (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 Figure 13. Star formation rate density at high metallicity (above solar, assuming Z / =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 andBisigello et al. (2018)-B18) and low mass slope of the redshift ∼ 0 MZR ( 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 . Note that starburst galaxies contribute to SFRD at relatively low metallicity and have little impact on the high metallicity cut of the SFRD. FMR based on the local MZR and SFMR. We introduce a parameter ∇ FMR0 to characterise 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 ∼ 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 prop-erties of the FMR on the parameters of the ∼ 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 (Z,z) models. We contrast the f (Z,z) derived that way with the f (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  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. budget). In addition to full f (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 3 than the redshiftdependent MZR extrapolated to high redshifts (as used in Chruslinska & Nelemans (2019)). However, while the high metallicity cut of the SFH at 3.5 is strongly affected by this assumption (Fig.  13), the low metallicity SFH estimate from  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 (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% solar metallicity at ∼0 ( ∼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 normalisation of the local MZR, while both MZR and starbursts strongly affect this quantity at > 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 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 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 < 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 ∼3. Model variations which account for the contribution of starbursts with B18/C17 based prescription lead to ∼2.5 times higher SFRD at 4 and shallower redshift evolution of the total SFRD at 2, more in line with the recent SFRD determinations at high redshifts. "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 will be made publicly available under this url: https:// ftp.science.ru.nl/astro/mchruslinska upon the acceptance of this manuscript. Figure A1. Dependence on parameter in eq. 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 / 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 / .

APPENDIX A: FMR VARIATIONS -ADDITIONAL FIGURES
In this Section we show additional figures illustrating the sensitivity of Z / (SFR,M * ) to the parameter ( Figure A1) and the Z / (SFR,M * ) dependence on various choices of the parameters of our observation-based model. Figure A2 shows the values of as a function of log 10 (SFR) obtained as described in Sec. 3.1.4 for different choices of the ∼ 0 MZR, ∇ FMR0 and SFMR used in this study. Figures A3 and A4 show the comparison between the Z / (SFR,M * ) obtained for different SFMR and ∇ FMR0 choices, respectively. 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. Figure B1 (B2) shows the fraction of the total SFRD that happens below 10% solar metallicity, i.e / < / -1 (above solar metallicity, i.e. > / ) as a function of redshift/lookback time, obtained for the different variations of our model discussed in Sec. 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.

APPENDIX C: Z / (SFR, M * ) IMPLEMENTATION WITH SFR-DEPENDENT ∇ FMR (AND )
As discussed in Sec. 3.1, Z / (SFR, M * ) at fixed SFR and M * < 0;SFR appear as lines roughly parallel to each other (see examples in, e.g. Mannucci et al. 2010, Telford et al. 2016, Cresci et al. 2019) and we assume that the Z / (SFR, M * ) slope is independent of SFR/M * . However, in their recent study Curti et al. (2020) fit the slopes of Z / (SFR, M * ) in individual log 10 (SFR) bins and indicate a possible dependence of on SFR. They find that decreases with decreasing SFR, reaching the value close to that of ∼0 MZR slope at low SFR (see figure 9 therein). When Z / (SFR, M * ) is extrapolated down to low (specific) SFR values, such dependence leads to considerably higher metallicity of galaxies with the lowest masses (M * 10 8 M ), and therefore can affect the low metallicity tail of the f SFR (Z,z) distribution obtained with our model. We explore additional variation of our model to illustrate the effect of such a dependence on our results. Within our description we can model a similar dependence of on SFR to that found by Curti et al. (2020) by implementing a SFR-dependent ∇ FMR (and therefore of , see e.q. 5) and requiring that the strength of the correlation weakens towards low SFRs, eventually disappearing at a certain low SFR (where approaches the slope of the z∼0 MZR). As discussed in Sec. 3.1.5, there is some observational evidence for the weakening of the SFR-metallicity correlation below a certain specific SFR threshold ( 10 −10 /yr, but it should be noted that this behaviour is difficult to constrain at low M * 10 8.5 M ). For simplicity, we assume that there is a linear dependence of ∇ FMR on log 10 (SFR), such that at log 10 (SFR) log 10 (SFR ∇0 ) ∇ FMR = ∇ FMR0 and at log 10 (SFR) log 10 (SFR no corr. ) ∇ FMR =0. We set log 10 (SFR no corr. )=-2, which roughly corresponds to the SFR of M * ∼ 10 8 M galaxies according to our z∼0 SFMR. We note that below that mass and SFR the observational FMR is not constrained 15 , and by assuming that for those galaxies the correlation disappears we obtain a conservative estimate of the low metallicity contribution of low mass galaxies. Results of our calculations are also sensitive to the value of the SFR threshold above which is fixed: the higher the value of log 10 (SFR ∇0 ), the lower the amount of the SFRD happening at low metallicity. The choice of this value is unclear and for the illustrative purposes we show our results using two example thresholds: log 10 (SFR ∇0 )=0 and log 10 (SFR ∇0 )=1. The resulting Z / (SFR, M * ) is illustrated in Fig. C1 using the latter threshold and for the different z∼0 MZR choices considered in our study. The overall f SFR (Z,z) distribution is only mildly affected by the explored variations with SFR-dependent ∇ FMR (see Fig. C2, the effect is analogous for the variations with the starburst implementation following Boco et al.). The redshift evolution of the peak metallicity is indistinguishable from that estimated with a fixed ∇ FMR = ∇ FMR0 case and the only visible difference is in the extent of the low metallicity tail of the distribution (the fixed ∇ FMR case leads to more extended low metallicity tail). With the low metallicity threshold of 10% solar metallicity (oxygen abundance) used in our study, the difference in the amount of the total SFRD happening at low metallicity estimated in the model variations with fixed and SFR-dependent ∇ FMR is amplified in the cases with flat z∼0 MZR (a MZR ∼0.3) and high MZR normalisation (compare lower edge of the dark brown range and the lower thick dark brown lines in Fig. C4 -note that the strength of this effect is very sensitive to the exact choice of the low metallicity threshold). In those cases the result is particularly sensitive to the assumed low mass/SFR FMR extrapolation, as at low redshifts only galaxies that are MZR outliers have metallicities falling below / < / -1 15 But note that galaxies with * 10 8 M and high specific SFR (i.e. SFMR outliers with high SFR) were found to follow the FMR (e.g. Mannucci et al. 2011;Hunt et al. 2016)   (see top left panel in Fig. C1). If there is no correlation between the SFR and metallicity at low SFRs, this leads to a negligible amount of SFRD happening below 10% solar metallicity at low to intermediate redshifts.
This paper has been typeset from a T E X/L A T E X file prepared by the author.     Table B3. Coefficients of the power law fits to the low metallicity cut of the cosmic SFH (lower and upper edges of the ranges as in Fig. 12) obtained for different variations of our observation-based model under the assumption of a non-evolving low mass end of the GSMF. Within each redshift range indicated in the first column total SFRD can be approximated with the dependence: SFRD=A· (1 + Figure C1. Same as Fig. 4, but assuming a SFR-dependent ∇ FMR (leading to SFR-dependent ). Top (bottom) panels assume S20 (PP04) ∼ 0 MZR. All panels assume SFMR with a moderate flattening at high masses (a =0.72 at the high mass end), log 10 (SFR ∇0 )=1 and ∇ FMR0 =0.27. Figure C2. Distribution of the cosmic SFRD at different metallicities and redshift for different Z / (SFR, M * ) implementations, shown for B18/C17 starbursts implementation (effect is analogous for the starburst implementation following Boco et al. (2021)). Brown contours and background colours correspond to Z / (SFR, M * ) with a fixed ∇ FMR0 = 0.27. Solid white (dashed black) contours correspond to Z / (SFR, M * ) modelled with ∇ FMR varying between log 10 (SFR)=-2 and log 10 (SFR)=1 (log 10 (SFR)=0) as described in Sec. C. Different 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 3 yr] (with the highest value corresponding to the innermost contour). The alternative Z / (SFR, M * ) implementation has a relatively small effect on the overall distribution, but leads to a less extended low metallicity tail. Figure C3. SFRD happening below 10% solar metallicity ( / < / -1) as a function of lookback time/redshift estimated with different Z / (SFR, M * ) implementations, shown for B18/C17 starbursts implementation (effect is analogous for the starburst implementation following Boco et al. (2021)). The light and dark brown coloured ranges are the same as in Fig. 12 and correspond to a fixed ∇ FMR0 = 0.27 and a MZR ≈0.3 and a MZR ≈0.6 respectively. Dotted gray/brown lines indicate upper and lower edges of the corresponding ranges obtained with ∇ FMR0 = 0.17. Solid (log 10 (SFR ∇0 )=0) and dashed (log 10 (SFR ∇0 )=1) lines indicate the upper and lower edges of those ranges obtained with a SFR-dependent ∇ FMR (leading to SFR-dependent , see Sec. C). Figure C4. Fraction of the total SFRD happening below 10% solar metallicity ( / < / -1) as a function of lookback time/redshift estimated with different Z / (SFR, M * ) implementations, shown for B18/C17 starbursts implementation (effect is analogous for the starburst implementation following Boco et al. (2021)). Top panel shows the estimates obtained for the low metallicity variations (i.e. corresponding to upper edges in Fig. C4), bottom panel shows the estimates obtained for the high metallicity variations. The green (purple) ranges span between the estimates obtained with Z / (SFR, M * ) with a fixed ∇ FMR0 = 0.3 and 0.17 and a MZR ≈0.3 (a MZR ≈0.6). Different lines indicate the corresponding estimate obtained with Z / (SFR, M * ) with a SFR-dependent ∇ FMR , varying between log 10 (SFR)=-2 and log 10 (SFR)=0 or 1 (see legend and Sec. C).