Metallicity of stars formed throughout the cosmic history based on the observational properties of star forming galaxies
Martyna Chruslinska Gijs Nelemans
astro-ph.GA


firstpage–lastpage
Calibration of quasi-static aberrations in exoplanet direct-imaging instruments with a Zernike phase-mask sensor
 III. On-sky validation in VLT/SPHERE
    
    A. Vigan1 
    M. N'Diaye2
    K. Dohlen1
    J.-F. Sauvage1,3
    J. Milli4
    G. Zins4
    C. Petit3
    
 
    Z. Wahhaj4
    F. Cantalloube5
    A. Caillat1
    A. Costille1
    J. Le Merrer1
    A. Carlotti6
    J.-L. Beuzit1
    D. Mouillet6

    Received 15 May 2019 / Accepted 18 July 2019
==================================================================================================================================================================================================================================================================





Metallicity is one of the crucial factors that determine stellar evolution.
To characterize the properties of stellar populations
one needs to know the fraction of stars forming at different metallicities.
Knowing how this fraction evolves over time is necessary e.g.
to estimate the rates of occurrence of any stellar evolution related phenomena 
(e.g. double compact object mergers, gamma ray bursts).
Such theoretical estimates can be confronted with observational limits to validate 
the assumptions about the evolution of the progenitor system leading to a certain transient.
However, to perform the comparison correctly one needs to know the uncertainties related to the assumed
star formation history and chemical evolution of the Universe. 


We combine the empirical scaling relations and other observational properties
of the star forming galaxies to construct the distribution of 
the cosmic star formation rate density at different metallicities and redshifts.
We address the question of uncertainty of this distribution due to currently unresolved questions, 
such as the absolute metallicity scale, the flattening in the star formation-mass relation or
the low mass end of the galaxy mass function.
We find that the fraction of stellar mass formed at metallicities <10% solar (>solar) since z=3 varies
by ∼18% (∼26%) between the extreme cases considered in our study.
This uncertainty stems primarily from the differences in the mass metallicity relations
obtained with different methods. 
We confront our results with the local core-collapse supernovae observations.
Our model is publicly available.



galaxies: abundances - galaxies: stellar content - galaxies: star formation -
stars: general - stars: abundances - stars: formation








§ INTRODUCTION

The distribution of the star formation over metallicities varies 
throughout the history of the Universe.
It is an important ingredient to estimate
the rate of occurrence of any stellar or binary evolution related phenomena,
such as different types of supernovae, double compact object mergers or gamma ray bursts
<cit.>.
Such theoretical estimates can be, for instance, confronted with observational determinations of rates
to validate the assumptions about the evolution of the progenitor star or system
leading to a certain type of transient. 
However, to perform the comparison correctly, one needs to know the uncertainties related to the
assumed star formation history and chemical evolution of the Universe. 
This is especially important in the case of transients whose
formation scenarios are particularly sensitive to metallicity, e.g. long gamma ray bursts (long GRB)
and stellar double black hole mergers, 
both forming much more efficiently at low metallicities 
(≲0.1 solar metallicity; see e.g. <cit.>, <cit.>, 
<cit.>, <cit.>
for long GRB and e.g. <cit.>, <cit.>, <cit.>,
<cit.>, <cit.> 
for double black hole mergers).


Recently, <cit.> 
demonstrated that the merger rate density
of double black holes estimated for the same description of the evolution of their progenitor systems
can be significantly different (even by a factor of ∼10) depending on the assumed
distribution of the cosmic star formation
rate density at different metallicities and time (probed by redshift; SFRD(Z,z)). 
This underlines the need for a better constrained picture
of the star formation history and chemical evolution of the Universe and of
understanding the associated uncertainties.


So far different groups have taken different approaches to determine the SFRD(Z,z),
often combining observations, theoretical inferences and/or cosmological simulations
 <cit.>.


In this paper we combine the empirical scaling relations from various observational studies 
describing the properties of star forming galaxies (star formation rate, mass, metallicity)
to construct the distribution of the cosmic star formation rate density at different
metallicities and time/redshift. 
Our method is outlined in Sec. <ref>.
We also address the question of the uncertainty of this distribution, given
the currently unresolved problems (discussed in detail in Sec. <ref>), 
such as the absolute metallicity scale <cit.> or the flattening (or lack of it)
at the high mass part of the star formation mass relation <cit.>.
Systematic evaluation of these uncertainties in the context of SFRD(Z,z) is still lacking in the literature.


Our results are summarized in Sec. <ref>.
In Sec. <ref> we apply them to calculate the volumetric core-collapse supernovae
rate and their local rate as a function of metallicity and contrast those quantities with observations.
We discuss the reliability of our results at high redshifts 
and present a brief comparison with simulations and earlier studies in Sec. <ref>.


The results of our calculations are publicly available at <https://ftp.science.ru.nl/astro/mchruslinska/>.
They can be applied to calculate the cosmological rates of various stellar evolution related
events and to asses their uncertainty due to uncertainties in observationally inferred SFRD(Z,z).
They can also be contrasted with the results from cosmological simulations.


Where appropriate we adopt a standard flat cosmology with Ω_M=0.3, Ω_Λ=0.7
and H_0 =70  km  s^-1 Mpc^-1 and assume a <cit.> initial mass function (IMF).



§ PROPERTIES OF STAR FORMING GALAXIES

While the chemical and star formation histories of individual galaxies are
heavily dependent on their environment and merger history, the picture
emerging from large galaxy surveys suggests that when a sufficiently large volume is considered,
the average properties of star forming galaxies follow relatively tight and simple, power-law like relations.
At a certain redshift both the star formation rate and metallicity correlate with the stellar masses
of galaxies <cit.>.
Those correlations span many orders of magnitude in mass and are present in the entire
redshift range where observations are available.


Below we provide an overview of the current observational results, 
underlining the open questions and uncertainties concerning the determination of various properties of star forming galaxies.
Based on the information presented in this section we decide on the relations and parameters used in our model.
Our choices are summarized below the relevant paragraphs.




 §.§  Galaxy stellar mass function of star forming galaxies 


The mass distribution of star forming galaxies can be inferred from the galaxy stellar mass function 
(GSMF, the number density of galaxies per logarithmic mass bin). 
Observationally such estimates have been provided by many groups
focusing on different redshift ranges:
e.g <cit.>, <cit.>, <cit.> at z∼0 ,
<cit.>, <cit.>, <cit.>,<cit.>,<cit.>,<cit.>
at intermediate and high redshifts.  
There are several studies aiming to determine 
the GSMF at redshifts as high z≳7 <cit.>
or even z≳8 <cit.>.


GSMF estimates for the same redshift bin can vary
substantially between different studies <cit.>, 
especially at the high and low mass end. 
Furthermore, the transition between redshift bins connecting results from different surveys
in general is not smooth, producing artificial jumps in the projected evolution of the GSMF over the cosmic history.
There is a number of factors that may be responsible for those differences.


For instance, stellar mass estimates usually rely on stellar population synthesis model fits to the
measured SEDs of galaxies. The result of SED-fitting depends on 
the parameters of the model (stellar population models, metallicity, dust law),
star formation histories and IMF used.
Different assumptions lead to systematic offsets in the mass estimates
<cit.>. 
Furthermore, GSMF estimates are based on different surveys, focusing on different redshift ranges (and wavelengths),
varying in depth and width of the covered sky area, completeness limits and introducing different biases.
On top of that, in case of GSMF for the star forming/quiescent subsample of galaxies
the final result is affected by the criteria used to select the interesting
population <cit.>. This choice is usually based on the colour - colour diagram, 
with various color indices and selection criteria used.


The observed (either for the active or total sample) galaxy stellar mass function generally declines with mass, 
shows a sharp cut-off at high masses (around M_∗∼10^10.7 at z∼0) 
and a power-law tail at low masses,
a relation well described and commonly fitted with a <cit.>(or sometimes double Schechter) 
function:

    n_gal=Φ(M_*)   d  M_* = Φ_* e^-M/M_*( M/M_cut  off)^α_GSMF dM_*

where n_gal is the number density of galaxies in a mass bin d M_*,
M_ cut off is the stellar mass at which the Schechter function bends, 
departing from a single power law with 
slope α_ GSMF at low masses to an exponential cut off at high masses. 
Φ_* provides the normalization (number density at M_ cut off).
The slope of the low mass end of GSMF is particularly weakly constrained. 
Even though the low mass galaxies are expected to be the most abundant in the Universe, 
they are also faint and difficult to observe especially at higher z. 
Also the mass completeness limit of any deep survey increases with redshift
(e.g. the sample from ZFOURGE survey used by <cit.> at z∼0.3 is 
complete down to M_∗∼10^8, 
while at z∼3 the completeness limit moves to M_∗≳ 10^9).
This leaves the GSMF in the low mass dwarf galaxy regime unconstrained.

Taking the fitted low mass slopes from different studies (see Tab. <ref>) at face value,
one would conclude that there is an overall tendency for the slope to steepen with z 
(α_ GSMF becomes more negative).
However, keeping in mind the differences in methods and surveys used by different authors,
correlations between the parameters of the fits 
<cit.>
and the fact that the source of discrepancies between the results from different studies for a single redshift
bin is not well understood, such a simple comparison may not reflect a true evolution of α_ GSMF.
Indeed, while some authors find evidence for α_ GSMF getting more negative with increasing z
<cit.>,
such evolution is not always found within one study covering a range of redshifts
<cit.>. 



Our choice: Instead of following the results of one group
we opt to average the estimates provided by different authors, 
similarly to the approach taken by <cit.> to constrain
their semi-analytic model.


The weakly constrained low mass end slope of the GSMF is treated separately.
We allow for two variations: one in which this slope is fixed to α_fix=-1.45 
as found by <cit.> (this seems to be a good compromise at least between z=0
and z=2 – during the bulk of the cosmic history, see Fig. <ref>) and one in which it increases
with redshift, following the linear fit shown in Fig. <ref>. 
See Sec. <ref> for the details.



 §.§ Metallicity
 

In this paper we use the word `metallicity` in a general sense 
(as a measure of the abundance of elements heavier than helium)
and use the symbol Z_O/H (or write explicitly 12 +  log (O/H ))
when referring specifically to the oxygen abundance ratio 
and Z when referring specifically to the mass fraction of heavy elements.

 The observations that we use provide estimates of metallicity in terms of oxygen abundance 
Z_O/H. 
However, the metal mass fraction is perhaps more useful for practical applications of our model
and hence we also present our results in terms of Z.
To convert Z_O/H to Z, we assume a simple scaling of the metal abundances with
Z_O/H that maintains the solar abundance ratios 
(i.e. log(Z/Z_⊙) = Z_O/H - Z_O/H ⊙).


However, there is little consensus in the literature regarding the value of
solar metallicity and solar composition.
Throughout the paper we assume solar abundances found by <cit.> 
 (Z_O/H ⊙=8.83 and Z_⊙=0.017), 
 since their results fall roughly in the middle of the range of the presently reported values
 (see appendix <ref>).
 All estimates of Z shown within our results were calculated assuming these values.
 We stress that this conversion is not unique and for the reasons discussed
 below in Sec. <ref> should be taken with caution.



  §.§.§ Oxygen vs iron abundance

The strong metallicity dependence of the efficiency of the
formation of various transients originating from massive stars 
(e.g. long gamma ray bursts, double black hole mergers)
is primarily driven by the abundance of iron.
This is because stellar winds from massive O-type and Wolf Rayet stars are driven 
by the radiation pressure on metal lines, 
and Fe easily dominates the atmospheric opacity due to its complex atomic structure
<cit.>.

Low wind mass loss rates are for instance necessary for a star to maintain high angular momentum,
which is needed to produce a long gamma ray burst <cit.>.
At the same time, lower wind mass loss allows for the formation of a more massive BH progenitor,
which then may form a BH in a direct collapse 
<cit.>.
The binary containing such a BH likely remains bound after its formation
and may evolve towards a merging double black hole system <cit.>.


Hence, it would be a more natural choice for our study to consider Fe instead of O abundances.
However, observational determination of the iron abundance is challenging and the number of
available results present in the literature much more limited than in the case of oxygen.
We thus rely on Z_O/H measurements assuming that it provides a good representation
of the overall metallicity (in particular the iron abundance) in the star forming material.


This is an important simplification to keep in mind, as the relative abundances of O and Fe in general do not
follow a simple proportionality, notably because 
the interstellar medium is enriched with different elements 
on different timescales <cit.>. 
While oxygen is generously released by massive stars (and hence on short timescales
∼10 Myr), a significant fraction of iron enrichment occurs via type Ia supernovae
(over much longer ∼Gyr timescales). 
As a result, young star forming systems often reveal an overabundance of oxygen relative to iron 
with respect to solar <cit.>.
The different timescales at which O and Fe abundance evolve in the interstellar medium
are reflected in the [Fe/O] vs [Fe/H] (or [O/H]) 
relation that can be obtained for a given stellar system <cit.>. 
The relation provides a better way to translate the O to Fe abundance for a particular system,
but is not universal (as it depends e.g. on the star formation history and the IMF) 
and hence is not applicable to our study.
Top-heavy IMF is another factor that can lead to increased ratio of 
α-elements to iron <cit.>.



 §.§ Mass – (gas) metallicity relation

The observed relationship between a galaxy's stellar mass and its metallicity
<cit.>
has been studied by many groups and using various methods to estimate metallicity.
The commonly employed set of methods make use of the optical emission lines coming from H II-regions
 to measure Z_O/H.


The most direct approach is to measure auroral lines, which allow to estimate the electron temperature
of the gas within the region, that in turn is a strong function of metallicity 
<cit.>.
However, auroral lines are typically weak,
 especially at high metallicities and cannot be used as tracers at Z_O/H≳ Z_O/H ⊙.
Furthermore, direct method based Z_O/H can be underestimated by even 0.4 dex in high-Z_O/H
environments if temperature gradients or fluctuations are present within the H II-region <cit.>.


To overcome those issues, several calibrations that allow to
translate ratios between the fluxes of strong emission lines into metallicity  have been
developed.
Those calibrations include empirical methods based on measurement of the electron temperature of the gas
<cit.>,
theoretical methods that rely on photoionization models 
<cit.>, 
or combination of the two <cit.>.
The calibration dependent methods (referred to as strong line methods) lead to large differences
between the measured Z_O/H, with the offsets of even ≲0.7 dex <cit.>.
The direct method and empirical calibrations typically lead to 2-3 times lower estimates than the theoretical ones, 
with combined methods falling in between. Those differences lead to different shapes and normalizations of the final MZR
<cit.>.

Several studies used metal recombination lines <cit.> that are weakly dependent on electron
temperature and insensitive to temperature fluctuations 
and may serve as an independent (next to emission line methods) indicator of metallicity.
However, observations of those faint lines are even more challenging than in the case of auroral lines.
Recombination line studies found metallicities consistent with those indicated by the 
calibrations based on photoionization models.


Alternatively, the metallicity of the star formation can be estimated with the help of young, massive stars.
At the evolutionary timescale needed to form supergiants (a few 10 Myr) their host galaxy's interstellar medium
experiences very little chemical evolution and hence those stars are perfect candidates for that purpose.
Spectra of both blue and red supergiants were used to measure metallicities in galaxies even beyond the Local Group
<cit.>.
Another approach to infer metallicity of the star formation using stars is to look at integrated light
spectra of young massive clusters with ages <100  Myr and typical masses ∼10^4
<cit.>, which are mostly found in highly star forming galaxies.
Both approaches lead to metallicity estimates consistent with the direct method measurements
at low metallicities and O3N2/N2 empirical calibrations as found by <cit.> at high Z_O/H, 
being ∼0.4 dex lower than the 
metallicity estimates from calibrations based on photoionization models
<cit.>.


The source of these discrepancies between different methods is presently not clear, 
which makes it very difficult to set an absolute metallicity.

 Our choice: 
Since the calibration leading to the correct estimation of the metallicity is not known,
we use the mass-metallicity relations based on 4 commonly used calibrations from the studies by:
<cit.> and refined by <cit.> (M09), <cit.> (T04), 
<cit.> (KK04) and <cit.> O3N2 calibration (PP04).

These relations cover the range of possible slopes and normalizations
for gas MZRs obtained with different methods to measure metallicity
<cit.>. See Sec. <ref> for the details.




  §.§.§ Evolution with redshift


In course of the history of the Universe next generations of stars form and evolve, 
gradually enriching the surrounding medium with metals. 
A fraction of this metal rich material can be lost from galaxies e.g. due to feedback 
from supernovae explosions or AGN activity or diluted by the inflowing metal-poor material.
Still, one would expect to see the imprint of this general enrichment as some form of 
the evolution in MZR with redshift.
The metallicity of distant galaxies can be estimated almost exclusively with
the strong emission line measurements.
Different emission line diagnostics are used at different redshifts, 
which may introduce artificial evolutionary trends or mask the true evolution of the relation.
Furthermore, those diagnostics are calibrated based on observations of local galaxies,
and it is not clear whether the conditions in which those lines are formed
do not change with z
<cit.>.
Reassuringly, <cit.> argue that these effects do not strongly affect the 
abundance estimates from nebular lines and <cit.> 
show that the commonly used diagnostics can be reliably applied up to z∼2.


Nonetheless, all of those issues make the study of the MZR evolution with redshift extremely difficult.
<cit.> study the MZR between z=0.05 and z=0.75 for a large sample of star forming galaxies
 with M_*>10^9.5. They use three different strong line theoretical calibrations
to estimate metallicities (McGaugh et al. (1991) et al.; Kobulnicky & Kewley (2004) and Tremonti et al. (2004)).
Within the considered mass and redshift range they find no evidence for a mass dependent evolution 
and see a clear indication of decrease in metallicity with z for all three calibrations.
However, the inferred rate of metallicity evolution is calibration dependent, with the slowest rate 
(- 0.16  dex/z) revealed when the <cit.> calibration is used, which also leads to the least steep MZR.
The other two calibrations lead to similar decrease rates (∼ - 0.26  dex/z).
<cit.> used the direct method to study the evolution of MZR at lower stellar masses
up to z∼1 and also concluded that at a fixed stellar mass
within that redshift range metallicity decreases by around 0.25 dex when compared with the 
z∼0 direct method based relation by <cit.>
and the shape of the relation at z≈0.5-1 is consistent with that found at z∼0.
On the other hand, the results obtained by <cit.> and <cit.> who studied metallicities
of active galaxies at z∼3 combined with z∼2 measurements from <cit.>
indicate that MZR evolves in a mass dependent way, at least in the high-mass part (M_∗≳ 10^9) 
probed in those studies. They find that the characteristic mass at which the MZR flattens increases with z 
and the overall evolution is stronger at lower masses.
A similar change in shape of the MZR was found by <cit.>. Those authors argue that the mass independent
evolution reported by <cit.> may be caused by the differences in sample selection and the
fact that their sample was limited to considerably higher mass galaxies than in <cit.>.


The MZR can be studied with the emission lines up to z∼3.5,
where the optical lines used to measure metallicity move out of the near IR bands
accessible to current ground-based spectrographs. 
Different tracers and methods need to be used beyond that redshift.
At cosmological distances metallicity can be estimated through absorption, 
e.g. using Damped Lyman Alpha systems <cit.>. 
However, it is not clear how the metallicity measured
that way relates to the metallicities of star forming galaxies and how the two methods should 
be compared <cit.>.
<cit.> used absorption lines in the interstellar medium of long gamma ray burst
host galaxies to infer the MZR up to z∼5 and concluded that the
relation continues to evolve at z>3, although it is not clear whether the long GRB host galaxies
provide unbiased sample of high redshift star forming galaxies.


Alternatively, one can use the information from the rest-frame UV spectra
that can be observed in the optical range for z≳ 2 
<cit.>.
<cit.> validated the correlation between the metallicity and the equivalent width of
absorption features in the rest-frame UV 
 at z∼2-3 and assuming that it also holds at higher z,
 applied that relation to probe the MZR at z∼5.
 They find a very weak correlation between the stellar mass and UV based metallicity at z∼5
within the probed mass/metallicity range, but their results come with large uncertainties
<cit.>.
If those findings can be directly confronted with the gas-phase MZR from <cit.>,
they show no evidence for evolution of the MZR since z∼3.5.
However, the authors stress that the uncertainties involved in the measurements of various parameters 
(in particular metallicity and stellar mass) need to be reduced to allow for any firm conclusions.

 Our choice: 
At z≤3.5 we interpolate between the literature results.
The evolution of MZR is unconstrained at higher z. 
We assume that the normalization continues to decrease, while the shape remains the same as at z∼3.5.
The rate of evolution depends on the metallicity calibration.

See Fig. <ref> and Sec. <ref> for the details.




  §.§.§ Scatter in the MZR
 

Metallicity as described by the MZR represents the average metal content of a galaxy
of a certain stellar mass and at a certain redshift.
<cit.> find that for a given stellar mass,
there is an intrinsic scatter of σ∼0.1 dex (containing 68% of the metallicity distribution)
around the best-fit mass - metallicity relation. <cit.> find a similar scatter
of 0.08 - 0.13 dex, depending on the metallicity calibration.

There is some evidence that the scatter may increase towards lower masses
<cit.>.
It is presently not clear whether the scatter in MZR evolves with redshift.
Determination of the intrinsic scatter is limited by the accuracy of the estimation
of observational uncertainties, which is more challenging in high redshift studies.
However, the analysis performed by <cit.> suggests that there is no significant
evolution in the magnitude of the MZR scatter up z∼0.8.

 Our choice: 
For a given stellar mass and redshift we assume that there is a normally distributed scatter around
the mean gas metallicity given by the MZR, with the dispersion σ_0=0.1 dex 
(the intrinsic scatter in the MZR). 
We allow σ_0 to increase linearly with decreasing mass at M_* < 10^9.5.




  §.§.§ Distribution of metallicity within galaxies


The average metallicity of a galaxy roughly corresponds to the metal content
that would be measured in this galaxy at the distance of ∼ r_e 
one effective radius from its center <cit.>. 
Due to metallicity gradients present within galaxies there is a certain
range of metallicities with which stars can form inside their host.
The azimuthal metallicity variations are typically negligible <cit.>.
However, typical radial abundance gradients in the local disk galaxies are around -0.1 dex/r_e 
(∼ -0.03 dex/kpc) <cit.>. 


The detailed studies of the distribution of metallicities of HII regions within
z∼0 disk galaxies carried out by <cit.>
and <cit.> with the CALIFA survey
show that the distribution of metallicities at which the star
formation proceeds is roughly symmetric with respect to the average
value found at ∼ r_e. 
However, the exact range of metallicities of the HII regions
and the metallicity gradient depend on the assumed
metallicity calibration <cit.>
[
For instance, 75% contours of the density distribution
of HII regions shown in Fig. 9 in <cit.> span the range of Δ_Z_O/H∼0.34 dex
in metallicity when the <cit.> O3N2 calibration is used.
<cit.> found a narrower range of Δ_Z_O/H∼0.22 dex.
However, the metallicity calibration used in this study <cit.> 
leads to shallower metallicity gradients.
In both cases the range is roughly symmetric with respect to the metallicity corresponding to the
value measured at r_e, with the HII regions in the inner parts of galaxies having metallicities
ranging from Z_O/H(r_e) to Z_O/H(r_e)+Δ_Z_O/H/2 and the outer parts from 
Z_O/H(r_e) to Z_O/H(r_e)-Δ_Z_O/H/2. 
].

 

Our choice: 
We assume that the distribution of metallicities at which the star formation proceeds
within galaxies can be represented by a normal distribution with dispersion σ_∇ Z_O/H=0.14 dex.
This particular choice of σ_∇ Z_O/H corresponds to half of the average range of metallicities
 of HII regions^2 found by <cit.> and <cit.>, 
 averaged between the two studies.

 See appendix <ref> for additional discussion.




 §.§ Star formation – mass relation


The stellar mass and star formation rate of star forming galaxies are strongly correlated,
giving rise to the star formation – mass relation 
<cit.> 
also called the star forming main sequence.


SFRs are measured based on the luminosities of galaxies
observed in certain bands that correlate with the recent (≲ 100 Myr)
star formation activity.
Those luminosities need to by corrected for dust attenuation, 
which is one of the main sources of uncertainty in SFR estimation.
Furthermore, conversion from luminosity to SFR is sensitive
to the assumed initial mass function (IMF) and metallicity.
The most common tracers of SFR are UV luminosity (usually at ∼1500 – 2800 Å),
certain recombination lines e.g. Hα and IR continuum emission between ∼ 3 - 1100 μm.
Different indicators are sensitive to the star formation on different timescales:
Hα line typically probes short ∼10 Myr timescales
while UV and IR provide information on SFR averaged over longer ∼100 Myr timescales
<cit.>.


The star formation main sequence at a certain redshift is commonly described as a power law relation
connecting the galaxy stellar mass M_∗ and SFR

    log(SFR/ yr^-1) = a × log(M_∗/) + b

with coefficients a and b describing slope and normalization of the relation respectively.
The typically found low-mass end slope values fall between 0.75 and 1 <cit.>.
Some studies find departure from the single power-law given by eq. <ref>,
identifying a flattening at the high mass end 
<cit.>,
while the other find no evidence for such a turnover <cit.>.
Hence, the presence and the amount of flattening in the SFMR is currently not clear.
<cit.> suggest that the presence of the turnover may depend on the method
used to filter out the quiescent galaxies from the sample used to construct the SFMR.

 Our choice:
As a base relation we adopt the recent result obtained by <cit.>,
who focused on the lower mass part of the SFMR where it can be described as a single power law.
Since the presence (and the degree) of the flattening at the high mass end of the relation is debatable,
we explore three variations of its shape: single power law at all masses (no flattening),
broken power law (moderate flattening) and almost constant SFR at high masses (sharp flattening).


The details can be found in Sec. <ref>.






  §.§.§ Evolution with redshift


<cit.> converted results from 25 studies probing M ≳ 10^10 and 
reaching up to redshift z∼5 to a common set of calibrations and
concluded that the high mass end slope of the SFMR  shows a mild
redshift evolution (see Fig. 8 therein). 
To our knowledge there is no similar analysis constraining the lower part of the SFMR.

Regardless of the precise form of the SFMR, its normalization is known to evolve with redshift.
This evolution is commonly parametrized by adding a factor 
c× log(1+z) to eq. <ref> with power law exponent c.
Values reported in the literature range from c∼1.8 up to c∼4 at redshifts z≲2
<cit.>.
Some studies find evidence for higher values of c at high mass part of the SFMR
<cit.> than at the low mass part.
Combining different datasets, <cit.> find c∼2.8.
At z≳2 the SFMR normalization shows little to no evolution
<cit.>
in tension with theoretical predictions <cit.>, based on which
one would expect the SFMR to decrease monotonically with cosmic time with c∼2.2-2.5.

 Our choice:
We assume that the low mass end slope of the SFMR does not evolve with redshift.
The evolution of the high mass part depends on the variation (see Sec. <ref>).
The normalization increases with redshift as c× log(1+z)
 with c=2.8 at z≤1.8 following <cit.>. 
 This value falls roughly in the middle of the range of values found in the literature.
 At higher redshifts we use c=1 to reproduce the observed flattening
 in the redshift evolution of SFMR normalization.
 The redshift at which we change the value of c corresponds to the redshift
 of the peak in the cosmic star formation history <cit.>.
 The redshift evolution of our SFMR is shown in Figure <ref>.




  §.§.§ Scatter in the relation

Similarly to the MZR, there is 'intrinsic' scatter (σ_ SFR) in the SFMR.
Its determination requires disentangling the
measurement error, redshift evolution within the sampled range
and the intrinsic scatter and proves challenging.
<cit.> found σ_ SFR=0.2 dex and
a similar value was recently obtained by <cit.>.
Those values are on the lower side of estimates present in the literature
which report the SFMR scatter about ∼0.3 - 0.4 dex 
<cit.>,
while <cit.> and <cit.> found σ_ SFR as high as 0.44 dex.
As argumented by <cit.>, those differences may be partially attributed to different SFR timescales
probed by different indicators.
They also point out that using the same data to derive SFR and stellar masses induces
correlation between the two which might artificially decrease the scatter found in the studies that do so.
<cit.> suggest that the scatter may increase towards lower stellar masses, 
but no clear trend neither with mass nor with redshift was found in other observational studies 
<cit.>.

 Our choice:
For any M_* and z we assume that the SFR is normally distributed around the mean given
by the SFMR with the dispersion σ_ SFR=0.3 dex.



 §.§ Fundamental metallicity relation


<cit.> suggested that there is a more general relation connecting
all three quantities characterizing galaxies discussed above: 
stellar mass, gas-phase metallicity and star formation rate, 
further investigated by <cit.> and called fundamental metallicity relation (FMR).
The relation is such that galaxies of the same stellar mass showing higher than average SFR 
also have lower metallicities.
A similar correlation was found in a number of studies
<cit.>.
The correlation weakens/ceases at very high stellar masses M_∗≳ 10^11.
<cit.> found no evolution in the FMR up to the redshift of z∼2.5, 
i.e galaxies at those redshifts
are found on the same M_*-metallicity-SFR 3D plane as the local galaxies, contrary to <cit.>
who found the redshift dependence of FMR.
The exact form of the observationally inferred FMR depends on the method used to select galaxies, 
measure metallicity, SFR and stellar masses.
Variations in α-enhancement may also influence the observational estimates of FMR <cit.>.

 Our choice:
The existence of the FMR, regardless of its precise functional form, means that
we cannot choose SFR and metallicity of a galaxy of a certain mass independently.
To account for the anti-correlation in SFR and metallicity for a given M_∗,
we assume the scatter in both relations is anti-correlated
[For instance, if a galaxy's SFR (chosen from a
normal distribution with σ_ SFR=0.3 dex and mean μ_ SFR given by the SFMR) 
is higher than indicated by the mean then the metallicity assigned to that galaxy will be lower
than the mean μ_ MZR resulting from the MZR
(and vice versa; e.g. if SFR(M_∗,z)=μ_ SFR+σ_ SFR/3 then
Z(M_∗,z)=μ_ MZR-σ_ 0/3).
]




 §.§ Initial stellar mass function


The IMF influences many of the observable properties of stellar populations and galaxies (e.g. SFR,M_*)
and hence also the empirical scaling relations.
In this paper we assume a non-evolving, universal IMF.
This is a common assumption and it is already introduced within the observational
studies whose results we use in our work.


However, the IMF has been theoretically predicted to vary with star-forming conditions
in such a way that in metal poor environments and/or regions with warmer gas the IMF becomes more top-heavy
<cit.>.
Studies based on the resolved stellar populations in Local Group galaxies
reveal no clear evidence for systematic trends in 
the IMF with metallicity (e.g. , but see ).
Some evidence for such variations has been found in studies of globular clusters
or ultra compact dwarf galaxies and within starburst galaxies <cit.>.
Furthermore, recent studies found radial IMF variations within massive early type galaxies, 
with the IMF generally found to be more bottom-heavy in their centers 
(e.g. , but see ).
While there is growing evidence for no general universality of the IMF,
it is presently not clear how the IMF varies with environment or redshift 
(and to what extent those variations can arise due to the differences in methodology).
Self-consistent modeling of the variable IMF is challenging
and the proposed models struggle to reproduce all observational constraint simultaneously
(e.g. , but see )


We note that significant variations of the IMF (e.g. with metallicity), if present,
could have a strong and non-straightforward effect on our results.
However, with the current state of knowledge a meaningful assessment of this effect 
is a complex (if not impossible) task and is beyond the scope of this study.

 

Our choice:
We assume a <cit.> IMF with the mass limits 0.1 - 100  (independent of the environment or redshift).
We correct for the differences between the IMFs adopted in various studies as described in the appendix
<ref>.




§ METHOD


To construct our model and calculate the distribution of star formation rate density across 
different metallicities and redshifts,
we combine the observationally inferred galaxy stellar mass function (GSMF) of star forming galaxies
with observational relations connecting stellar masses of galaxies with their gas-phase metallicity (Z_O/H)
and star formation rate (SFR). Our method is schematically summarized in Figure <ref>.

The procedure can be summarized as follows:

 
  * We consider galaxies with stellar masses between 10^6 - 10^12
  and divide this mass range into N_gal mass bins Δ M_j equally spaced in logarithm.
  Within each redshift/time bin we calculate the number density of galaxies falling within
  each mass bin n_Δ M;j, using the galaxy stellar mass function of star forming galaxies.
  This mass density is assumed to be constant within the redshift/time bin.
 
  * From each mass bin we randomly choose N_sampl stellar masses. 
  Each of those masses M_i is associated with a certain number density
  of galaxies n_M_i=n_Δ M;j/N_sampl
 
  * Each mass M_i is then assigned a certain star formation rate SFR_i, 
  drawn from a SFMR for a given redshift and taking into account the scatter in the relation. 
  This SFR is assumed to be constant within the redshift/time bin.
 
  * For each mass M_i we draw the corresponding gas metallicity Z_O/H  i from the MZR 
  at a given redshift, taking into account
  the scatter in the relation and correlation between the SFR and metallicity.
   Again, Z_O/H  i is assumed to be constant within the redshift/time bin.
  
  * We then calculate the amount of mass formed in stars within the 
  time/redshift bin per unit of comoving volume in galaxies of different masses
   and metallicities
 and sum the contributions from galaxies that fall within the same metallicity bin.

We use different variations of the assumptions to estimate the uncertainties of our results.
In the following subsections we provide the details of the construction of our model and
describe how the different observational results introduced in Sec. <ref> are combined.







 §.§  Galaxy stellar mass function
	

We use the published best-fit parameters to the Schechter (or double Schechter) function describing
GSMF of star forming galaxies in different redshift bins at z≲6. 
The studies used in the analysis are listed in Table <ref>.
The last column of Table <ref> gives the limits of the redshift bins 
in which the GSMF was estimated in each of those studies.
We choose several values of redshift z_ avg
that overlap with the redshift bins in at least three of those studies
[
with the exception of z=5 where we use only one study by <cit.>, 
as this is the only result constraining GSMF of active galaxies at z>4.
]. 
Those values are listed in the first column of Table <ref>.


To obtain the number density of star forming galaxies of a certain mass
and at a certain redshift we:


(i) calculate the number density for that mass at z_ avg redshifts
as the average from the number densities calculated using the Schechter function fits
referenced in Tab. <ref>


(ii) interpolate between the values found at different redshifts z_ avg


We treat the low and high mass part of the GSMF separately, 
calculating the number density of galaxies using the full Schechter fits as
described above at M_∗>M_ GSMF and assuming that at lower masses
the GSMF of star forming galaxies is described as a single power law with slope α_ fix.
We allow for two variations: 

 
  * α_ fix=-1.45 at all redshifts
 
  * α_ fix=α_ fix(z) steepening with 
      redshift according to the linear fit shown in Fig. <ref> 
      up to z=8 and fix α_ fix=-2.14 at z>8 [
      The high redshift study by <cit.> shows little to no evolution of the GSMF at z≳8
      ].

We assume that the mass separating the low and high mass part of the GSMF increases with redshift as
log( M_ GSMF/) = 7.8 + 0.4 z
[
This roughly corresponds to the mass completeness limit as a function of redshift from <cit.>
(see fig. 2 therein). This limit in other z<5 studies is typically higher, thus we allow for some
extrapolation of the fitted Schechter function beyond that limit before fixing the slope.
]
up to z=5 and equals log( M_ GSMF/)=9.8 at
higher redshifts.


At z<0.05 we use the estimate obtained for z=0.05.
To our knowledge there are no measurements of GSMF for the star forming population of galaxies at z>6.
To extend our analysis to higher redshifts and constrain the GSMF evolution, 
we include the recent results obtained for the total sample of galaxies
at z∼7, z∼8 and z∼9 <cit.>. 
This is justified as at such high redshifts the fraction of passive galaxies is expected to be very low
<cit.>.
At z>9 we assume that the normalization of the GSMF continues to mildly decrease
at the rate found between z=8 and z=9.
The resulting GSMF is shown in Fig. <ref>.





 
 



 §.§ Mass-(gas) metallicity relation
 

We start from the mass metallicity relation fitted by <cit.> 
(in 3 redshift bins z∼ 0.07, 0.7 and 2.2)
and refined by <cit.> (at z∼3.5)
[
We convert the results from <cit.> to <cit.> IMF used by <cit.> 
by applying the correction suggested by these authors and then increase the
logarithm of the mass by 0.03 dex to convert these results to <cit.> IMF.
]
We use a different parametrization of the MZR than in those two studies, given by eq. <ref>
and proposed by <cit.>.

    12 + log[O/H] = Z_O/H  asym - log[ 1 + (M_∗/M_TO)^-γ]

This parametrization allows to avoid an artificial turn-off at the high mass part of the relation.
The parameters are: γ describing the low–mass end slope,
M_ TO – the mass at which the relation begins to turn/flatten and Z_O/H  asym – the
asymptotic metallicity of the high-mass end.
The relation given by <cit.> was refitted in each of the redshift bins using this 
parametrization.
We refer to this form of MZR as M09 throughout the paper.
The gas metallicities in <cit.> at 12+log[O/H]>8.35 dex were obtained with the
 Kewley & Dopita (2002) strong line calibration.
We translate the base M09 relation at each redshift
to three other calibrations 
(<cit.> – T04, <cit.> – KK04 and <cit.> O3N2 – PP04) 
applying the conversion method described in <cit.> at the 
high metallicity part (12+log[O/H]>8.35 dex) and fit the relation given by eq. <ref>
to the converted MZRs.

The fitted parameters are summarized in Table <ref>
[The slope of the PP04 MZR was fixed at z>2.2 to the value fitted for that redshift 
(could not be fitted with M09 data at metallicities higher than 8.35 alone)].


We interpolate between the redshift bins to obtain the MZR at 0 ≤ z ≤ 3.5 for each calibration.
We assume that at higher redshifts the relation evolves in a mass independent fashion, i.e. only the
normalization decreases. The rate of decrease in normalization (dZ_O/H/dz)
at z>3.5 is assumed to be the same as the decrease in metallicity at M_∗=10^11 between 
the two highest redshift bins (z=2.2 and z=3.5).
Hence, for a given calibration the MZR at z>3.5 is given by:

    (12 + log[O/H])_z>3.5 = (12 + log[O/H])_z=3.5 + dZ_O/H/dz × (z-3.5)

For a galaxy with mass M_* at redshift z, we draw its average metallicity
from the normal distribution centered at the metallicity taken from the MZR
and with dispersion σ_0=0.1 dex at M_*>10^9.5 
and σ_0=–0.04 log(M*/)+0.48 dex at higher masses.

We add a normally distributed scatter around this value with σ_∇ Z_O/H=0.14 dex
to account for the distribution of metallicities at which the stars are forming within their host galaxies.




 

  


 §.§ Star formation - mass relation
 
 Following <cit.>, we describe the low mass part of the SFMR 
 as a single slope power law with a=0.83 (see eq. <ref>).
 This slope does not evolve with redshift.
 We explore three variations of the SFMR shape at high masses:

 
  
  * no flattening - we use the relation given by <cit.> 
  in the entire mass range (extrapolating beyond the mass range
  7<log(M_∗)<10.5 covered in their study when necessary)
  
  * moderate flattening - we use the relation given by <cit.>
	at low masses and modify the slope at M_∗≥ 10^9.7
	according to the results obtained by <cit.> for  M_∗≥ 10^9.7
	and 0<z<6. 
	In this variation the slope 
	of the high mass part of the SFMR is less steep than that of the low mass part and
	evolves with redshift, becoming steeper with increasing z.
	At z>6 we extrapolate the high mass-end slope evolution with time found by these authors.
  
  * sharp flattening - we use the relation given by <cit.>
	at low masses and combine it with the relation obtained by <cit.>
	[
	<cit.> used a different parametrization of the SFMR (as opposed to the one given 
	by eq. <ref>):
	log( SFR ) = s_0 - log[ 1+ ( M_∗/ M_TO; SFMR)^-γ]
	,where γ defines the low mass-end slope, M_TO; SFMR is the stellar mass above which the
	relation begins to flatten and asymptotically approaches a peak value s_0
	]
	, which leads to flattening in SFMR at stellar masses above a 
	certain turnover mass M_∗≥ M_TO; SFMR.
	In this case SFR is almost constant with increasing M_∗ at M_∗≥ M_TO; SFMR,
	as opposed to the `moderate flattening` variation, 
	in which the slope changes but the SFR is still increasing with increasing mass.
	<cit.> concluded that M_TO; SFMR increases with redshift in the redshift range
	0.5<z<4 covered in their study. At z<0.5 (z>3.28
	[
	We use the value of M_TO; SFMR at z=3.28 instead of that corresponding to the edge of 
	the last redshift bin (3<z<4) used in the study by <cit.>, because their polynomial fit 
	log(M_TO; SFMR)=0.9458 + 0.865z - 0.132z^2 reaches maximum at z≈3.28 
	and a further decrease in M_TO; SFMR seems to be an artifact.
	]) 
	we fix M_TO; SFMR to the value that results from
	the relation fitted by these authors at z=0.5 (z=3.28).   
 
 
 
 In all variations we keep the normalization given by <cit.> 
 (converted to <cit.> IMF) at low masses and z∼0
 and tie the high mass part of the relation to the low mass part accordingly.
 The normalization increases with redshift as c× log(1+z)
 with c=2.8 at z≤1.8 and c=1 at z>1.8 (see Fig. <ref>).
 

 Similarly to the MZR, for any given stellar mass and redshift we 
 assume that the SFR is normally distributed around the mean given by the SFMR
 with σ_ SFR=0.3 dex.



 §.§ The model and its variations


Table <ref> summarizes the default choice of parameters used in our calculations.
The parameters that were varied are indicated in bold. We perform the calculations for 24 variations
of the parameters of the model. Each variation is defined by the choice of MZR, SFMR flattening and α_fix.

 



§ RESULTS


One of the goals of this study is to find the observation-based 
distribution of the cosmic star formation rate density at different
metallicities and redshift SFRD(Z,z). 
Given the currently unresolved observational issues, 
especially the unknown source of the differences 
between the metallicities measured using different calibrations
and the shape of the high mass part of the SFMR (see sec <ref>), 
there is no simple, single answer to the question about the shape of the SFRD(Z,z).
Hence, instead of indicating one `best` SFRD(Z,z) we first choose an example `moderate`
variation of our model to discuss the general characteristics of the resulting distribution. 
We then discuss the differences between the variations allowed by the current observations and 
focus on the extreme cases that lead to the highest fraction of 
low/high metallicity star formation,
which delineate the uncertainty of our result.



 §.§ The distribution of the cosmic star formation rate at different
 metallicities and redshift: general picture


 



Figure <ref> shows distribution of the star formation rate density 
at different metallicities and redshifts calculated assuming 
a moderate flattening at the high mass end of SFMR, 
M09 mass metallicity relation and fixed low mass end slope
of the galaxy stellar mass function α_fix=-1.45.
We use this variation as an example to point out the main characteristics of the SFRD(Z,z)
that are common to all model variations:


 
  * At low z the star formation is concentrated at relatively high metallicities;
 the location of the peak (thick brown line in Fig. <ref>) depends primarily on the MZR

  * The metallicity at which SFRD(z)[SFRD(Z,z) integrated over metallicities]
 peaks decreases towards higher z, reflecting the evolution of the MZR.
 The decrease rate increases around the redshift ≳2 
 [At z≳2 the peak decreases almost linearly, 
 at a slightly higher rate than the normalization of the MZR
(-0.28 dex versus dZ_O/H/dz = -0.26 dex in the case shown in Fig. <ref>)]
   and depends on the MZR.

  * The global peak of the distribution is reached at z∼1.8, 
      coinciding with the peak of the star formation history of the Universe,
      and at Z that depends on the model variation (primarily MZR)

  * The distribution is asymmetric with respect to the peak metallicity,
      showing the extended tail at low Z (see also Fig. <ref>).
      The contribution from the tail to the total SFRD(z) at high z 
      is noticeably higher in the variations with α_fix=α_fix(z)
      (see Fig. <ref>).

The bottom panel of Figure <ref> shows the redshift evolution of the
ratio of the total SFRD(z) in our simulations 
to the SFRD from observational studies by
<cit.>, <cit.> an <cit.>.
We note that the only way the total observed SFRD was used in our model is to
set the redshift at which we change the rate of redshift evolution of the SFMR normalization
[
The observational studies discussed in sec. <ref> 
show that this transition occurs somewhere between z∼1-2
and here we choose z=1.8, which is an average redshift of the peak of the star formation history
of the Universe resulting from studies by <cit.>, <cit.> an <cit.>
]
.
Thus, this comparison is still a valid test for our model.
The total simulated SFRD shows a remarkable agreement with observations up to very high redshifts,
staying within a factor of 2-3 from the observed SFRD up to z∼8 for all variations that assume
α_fix=-1.45 (but not for the cases with α_fix=α_fix(z), see Sec. <ref>).

 Note that the model variations that differ in the assumptions about the
SFMR and/or GSMF lead to different total SFRD(z) (which sets the normalization of the SFRD(Z,z)).
In general, SFRD increases up to z∼ 1.8, decreases towards higher z and declines sharply around z∼ 8.
This evolution is a result of the interplay between the evolving SFMR and GSMF,
mostly their normalizations.
The shape of both relations also changes with redshift, 
but it has a secondary effect on the evolution of the total SFRD 
(except for the cases with α_fix=α_fix(z) at z≳ 4).
At z∼8 the normalization of the GSMF decreases abruptly (see fig. <ref>),
which results in an apparent jump in the normalization in Figure <ref>.




  §.§.§ The variations with α_fix=α_fix(z)




When the low mass end slope of the GSMF is allowed to evolve (steepen) with redshift according 
to the fit shown in Fig. <ref>, the resulting SFRD(Z,z) distribution starts to
deviate significantly from the picture shown in Fig. <ref> at z≳ 4.
The SFRD(Z,z) distribution calculated assuming the same relations as in the case
shown in Fig. <ref>, but with evolving α_fix=α_fix(z)
is shown in Figure <ref> (the other model variations with α_fix=α_fix(z) 
show the same qualitative behavior).


The normalization of the SFMR continuously increases with redshift, 
which means that galaxies at all masses produce stars at higher rates. 
At the same time, the normalization of the GSMF decreases and the cut off mass shifts
to smaller values, both reducing the number density of the most massive galaxies 
(with the highest SFR and Z at a certain redshift). 
α_fix decreases (the slope steepens), increasing the number of low mass galaxies
and the shape of the GSMF evolves almost into a single power-law (see fig. <ref>).
This, together with continuously increasing SFR(M_∗) produces a plateau in the total SFRD(z)
in our simulations at z∼4-7 (see e.g. fig. <ref>,
 showing the core collapse supernovae rate, which is proportional to the total SFRD(z)). 
The low mass galaxies begin to dominate the total star formation budget and the peak of
the SFRD(Z,z) shifts to very low metallicities (below 1%  at z∼5).


Such plateau in SFRD(z) at z∼4-7 is not found within the current observational studies 
(predicting a continuous decrease in the total SFRD(z) at these high redshifts)
and hence the bottom panel of Fig. <ref> shows that our simulations
start to significantly overpredict the total SFRD at z∼6. 


Taking this into account, we consider the high redshift (z≳ 4) 
predictions of the model variations with α_fix=α_fix(z) unreliable
and exclude them from further analysis, while still taking them into consideration at lower redshifts.
We return to that discrepancy at high redshifts in Sec. <ref>.



 §.§ The extreme cases
 

The main purpose of this study is to identify the model variations 
 that lead to the most extreme
pictures of the SFRD(Z,z) and hence delineate its uncertainty.

Different assumptions about the amount of star formation occurring
at low metallicities can lead to significant differences e.g. in the obtained properties of
the progenitors of merging double black holes and long gamma ray bursts and the inferred rates
of transients connected to those objects.
Keeping in mind the importance of metallicity,
we interpret extremes as variations that lead to the 
smallest and the highest fraction of stellar mass forming at low/high metallicity.

The main differences in the SFRD(Z,z) distributions obtained for different
model variations are the following:

 
  * The location of the peak metallicity and the slope of the curve indicating
       the peak metallicity at each redshift
 
  * The extent of the low metallicity tail and its contribution to the total SFRD(z)
 
  * Normalization of the distribution / total SFRD(z) 

Focusing on the fraction of stellar mass formed at low and high metallicity
means that we set aside the differences in normalizations (but see Sec. <ref>).

 
The location of the peak metallicity and the rate of its decrease with z are set primarily by the choice of the MZR
– the relation with high, slowly decreasing normalization will maximize the fraction of stellar mass forming at high metallicity. 
The slope of the MZR affects the low metallicity tail of the SFRD(Z,z) - the flatter the relation,
the smaller the difference in metallicity of stars forming in galaxies of different masses
and the the low-Z tail is reduced.


Furthermore, the choice of the SFMR regulates the contribution of the most massive galaxies (forming stars at high metallicities)
to the total SFRD(z). Hence, the SFMR with no flattening maximizes the high metallicity star formation,
and the reverse is true for the SFMR with sharp flattening


Finally, the low mass end slope of the GSMF regulates the contribution of galaxies with the lowest masses.
The steeper the slope, the more low mass galaxies 
(and hence more stars forming in the low metallicity tail of the distribution).


Thus, one can expect the extreme variations to be represented by the combination of
KK04 MZR, SFMR with no flattening and α_fix=-1.45 (the high metallicity extreme)
and PP04 MZR, SFMR with sharp flattening
[
The choice of the α_fix variation in the low metallicity extreme is not obvious.
However, as can be seen by comparing Fig. <ref> and <ref>,
the choice of the α_fix variation affects the SFRD(Z,z) mostly at z≳3, up to which
redshift both variations give similar results.
At the same time, the variations with α_fix(z) are excluded at z≳4 (Sec. <ref>)
and we focus on the discussion of the variations with α_fix=-1.45. See Sec. <ref>.
]
(the low metallicity extreme).



  §.§.§ SFRD(Z,z) for the extreme cases




The distributions of the star formation rate density across metallicities
and redshifts for the extreme cases identified in Sec. <ref>
are shown in Figure <ref>. 
The difference between the two panels visually demonstrates the uncertainty of the SFRD(Z,z)
due to unresolved questions in the determination of various characteristics of galaxies 
(absolute metallicity scale, high mass end of the SFMR and low mass end of the GSMF).

 To facilitate the comparison, we overplot the constant SFRD contours and the line
showing peak metallicity at each redshift for the moderate variation shown in Fig. <ref>.
The left panel of Figure <ref> shows the high metallicity extreme.
The low metallicity extreme is shown in the right panel.


As shown in the right panel in Fig. <ref>, 
except for the ∼0.4 dex offset, the redshift evolution of 
the peak metallicity for both extremes is very similar.
This stems from the fact that MZRs obtained for different metallicity calibrations
have similar shape (i.e. are relatively flat) in the high mass part 
that contributes the most to the total star formation.
At higher redshifts z≳ 2 the peak metallicity lines for both cases start to deviate
due to different rates at which the MZR normalization decreases for both relations
(the PP04 MZR normalization decreases at a higher rate than KK04; see Table <ref>).
The peak metallicity reaches our conventional low metallicity limit of 0.1
at z=5.8 and z=3.6 for the high and low metallicity extreme respectively
and never drops below  0.01 in the former case.


Note that the low metallicity extreme 
has a lower SFRD(z) and leads to smaller total stellar mass than the moderate variation
and the high metallicity extreme (see Sec. <ref>).




  §.§.§ Quantifying the differences between the variations

In this section we aim to quantify the differences between the model variations discussed above,
in particular to find the maximum difference between the fraction of stellar
mass formed at low and high metallicity across the considered variations.


For presentation purposes we adopt the conventional limits of 
Z<0.1 and Z> to define the low and high metallicity respectively
 <cit.>.
 We note that this particular choice does not affect our conclusions.
 

 In our analysis we focus on three redshifts ranges: z<0.5 to which we refer as the `local Universe`
and which can be probed with the current network of ground-based gravitational wave detectors, 
z<3 up to which redshift our model is still reasonably well backed up by the current observational
results and which captures the great majority of the star formation history of the Universe
and z<10 where we extend our calculations beyond redshifts directly probed in observational
studies relying on extrapolations.
Those redshifts correspond to ∼12.8 Gyr,  ∼2.1 Gyr and ∼0.46 Gyr 
respectively in terms of the age of the Universe.
The example distribution of stellar mass formed since z=3 across metallicities
is shown in appendix <ref>.

 Figure <ref> shows the fraction of stellar mass formed at low metallicities 
versus that formed at high metallicities since z=0.5, 3 and 10 for our model variations.
It can be seen that, as expected, variations leading to the highest fraction of mass formed at low metallicities
at the same time produce the lowest fraction of high metallicity stars (and vice versa) at all redshifts
and the extreme cases introduced earlier in Sec. <ref> can be easily identified
(low-Z extreme - brown circles; high-Z extreme - red squares).
Note that when the stellar mass formed since z>1.5 is considered, 
the fraction of stellar mass formed at low-Z is additionally increased 
if the low mass end slope of the GSMF is allowed to steepen with z
[
At redshifts z<1.1 α_fix(z) resulting from the fit shown in Fig. <ref>
is bigger than α_fix=-1.45 that we use in the variations in which the slope does not evolve with z.
Hence, there are less low mass galaxies in the variations with α_fix(z) in that redshift range and
a smaller fraction of stellar mass forms at low metallicities. For the variation with 
PP04 MZR and sharp flattening in the SFMR those fractions equalize at z∼1.5.
] (compare brown circles with and without edge).
The impact of different assumptions on the low-Z mass fraction is further discussed in appendix <ref>.

 The difference between the fraction of the mass formed at low and high metallicity
produced in the two extremes is quite significant.
Those values (summarized in Table <ref>)
range between ∼3%-15% and ∼4%-55% at z≤0.5,
∼9%-27% and  ∼1%-27% at z≤3 and
∼17%-43% and ∼0.7%-19% at z≤10 for the low and high metallicity fractions respectively.



For the example variation shown in Fig. <ref>
(with M09 MZR, moderate SFMR flattening and non-evolving α_fix)
the corresponding values are 6% and 38% at z≤0.5,
15% and 15% at z≤3, 
and  29% and 11% at z≤10.
We note that it does not fall exactly in between the low and
high metallicity star formation extremes,
but produces a higher fraction of mass at high metallicities with respect to a simple average.
While it represents the `most moderate` case considered in this study 
(and hence we refer to it as the moderate variation), 
the results obtained for the T04 MZR (for the same assumptions about the SFMR and GSMF) are comparable.
This highlights the question of degeneracies involved in our calculations.
Similar results could be obtained using another MZR (resulting from a different metallicity calibration),
or by combining it with different assumptions about the SFMR and/or GSMF 
that fall between the extreme cases considered in our work.

 In nearly all cases, the summed fraction of stellar mass produced at low and high metallicity is smaller than
0.5. This means that in those variations most of the stellar mass forms outside the assumed low and 
high metallicity regimes
(i.e. between 0.1 and  or, equivalently, between Z_O/H=7.83 and Z_O/H=8.83)
[
If only the local redshift range is considered, the exceptions are the variations with KK04 and
no or moderate SFMR flattening and T04 MZR with no SFMR flattening. 
We find that at z<0.5 in those cases more than half of the stellar mass
formed at metallicities higher than the solar value].
We look in more detail at the location of the peak of the distribution
- i.e. metallicity at which most of the stellar mass formed in different redshift ranges
and model variations in appendix <ref>.








 §.§ Total stellar mass vs the mass fraction


 

The different variations do not only lead to different fractions of stellar mass formed at high and low metallicity, 
but also produce different total stellar mass (see Table <ref>).
Variations with no flattening in the SFMR lead the highest total M_* because they
produce more stars in the high-mass galaxies (with high SFR and Z).


This effect reduces the difference in the total mass formed at low metallicity 
and increases the difference in terms of 
the total mass formed at high metallicity between the variations with different assumptions about the SFMR.


Using the values given in Table <ref>, we find that
the total mass formed at low metallicity across the model variations differs by less than a factor of 3
and is relatively well constrained within our study.

The corresponding values of the mass formed at high metallicity 
span a much wider range and the differences between the model variations reach up to a factor of ∼40.




§ APPLICATION OF THE MODEL: CORE COLLAPSE SUPERNOVAE RATES




Core-collapse supernovae (CCSN; i.e. Type II and Ib/c SNe) originate from short-lived massive stars
and can serve as an independent tracer of the star formation <cit.>.
We can thus confront the predictions of our model with the volumetric CCSN rate estimates at different redshifts,
assuming that the volumetric CCSN rate is proportional to the SFRD(z).
However, one should keep in mind that there is a number of factors that make this comparison not straightforward
and may cause deviations from the simple proportionality (see appendix <ref>).


In Figure <ref> we plot the observational CCSN rate estimates at different redshifts
against the star formation rate density for the two extreme cases and the moderate variation of our model
discussed in Sec. <ref>, multiplied by a constant factor k_CCSN∼ 0.009  ^-1.
This factor comes from a simple IMF-based scaling and corresponds to the 
efficiency of formation of CCSN progenitors assuming our default <cit.> IMF and CCSN forming
form single stars with masses between 8-25  
[
The k_CCSN for z≲1 and the moderate variation of our model
lies between k_CCSN∼0.0082 - 0.0097 ^-1
(depending on the exact redshift range and data sample considered in the fitting).
Assuming a fixed low mass end of 8  and our default <cit.> IMF, 
this corresponds to the upper CCSN progenitor mass end between ∼20 - 30 .
]
, which given the large scatter in the observations,
provides a satisfying match to the data.

 The knowledge of the distribution of the cosmic star formation over metallicities
allows us to compare not only with the global volumetric CCSN rate discussed above,
but also to perform a more detailed comparison with the CCSN rates measured as a function of metallicity.
We calculate the specific CCSN rate as a function of metallicity
(number of CCSN in a certain metallicity bin divided by the stellar mass
in galaxies in that metallicity bin)
and use the mass range for the CCSN progenitors estimated based on the comparison with the volumetric rates
to get the normalization.
Those rates can be confronted with the specific CCSN rates based on the Lick Observatory Supernova Search 
(LOSS) sample from <cit.>
[
The CCSN specific rates were obtained by summing the specific rates for type II and stripped envelope SN.
The rates given in <cit.> were estimated using different metallicity bins for the two types of SN. 
Here we use three metallicity bins (in the Tremonti et al. (2004) scale): 
8.71^+0.18_-0.24, 9.07^+0.04_-0.06, 9.17^+0.07_-0.04 
in which the corresponding 
type II SN rates are: 3.47^+1.06_-0.83, 0.90^+0.28_-0.22, 0.52^+0.15_-0.12
and stripped envelope SN rates are: 0.97^+0.66_-0.42, 0.34^+0.27_-0.16, 0.25^+0.17_-0.11
(Graur, private communication)
].
Figure <ref> shows the specific rates at z≲0.02, which roughly corresponds
to the redshift range relevant for the LOSS SN sample. We show the results for the model variations 
using the T04 MZR, as this was the metallicity calibration used in <cit.>.
[
We note that the choice of the MZR does not affect the result of the comparison.
To compare the observations with the results obtained for the model variation with a different MZR,
one has to convert the data points to a different metallicity calibration.
Applying the conversion method described in <cit.> to the observations
results in a shift in the metallicity of the data points which is consistent with the shift of
the model curves obtained for the different MZR and hence the conclusions drawn from the comparison do not change.
]


The metallicity measurements used in <cit.> are based on the SDSS observations and 
correspond to the central metallicity within the CCSN host galaxy.
The central metallicity is typically higher than the metallicity at the CCSN location 
(due to radial metallicity gradients, see Sec. <ref>).
It is also higher than the metallicity probed by the MZR, which gives the `averaged` global value
(and corresponds roughly to the metallicity that would be measured at the distance of one effective radius r_e).
Therefore, an offset in metallicity between the data and the results from our model can be expected.
To be conservative in our comparison, we extend the horizontal errorbars of the data 
shown in Figure <ref> by 0.24 dex towards lower metallicity,
assuming that the great majority of the observed CCSN explosion site metallicities would 
be contained within that range (see appendix <ref>).
This wide range, together with the uncertain translation of the
amount of star formation to the number of core collapse supernovae (which affects the location
of the model curves on the vertical axis) precludes us from drawing any strong conclusions from that comparison.
In general, we find a reasonably good agreement of our model with the data.
However, we note that the dotted curve in Fig. <ref> 
which corresponds to the model variation with the sharp flattening in the high mass end of the SFMR, 
falls slightly outside the extended uncertainty regions of the data points.
Hence, the model variations with a sharp flattening in the SFMR
may underestimate the amount high metallicity star formation in the local Universe.


We note that the slope of the high metallicity part of the relation plotted in 
Fig. <ref>
depends on the combination of the scatter present in the MZR (σ_0) and on the 
width of the distribution of metallicities at which the stars form within galaxies 
(σ_∇ O/H; see appendix <ref>).
The slope cannot be directly compared with the data from <cit.>, 
as they only provide the central metallicities of CCSN host galaxies,
but could potentially be probed with similar observations that provide metallicity measurements at SN location.







§ DISCUSSION

In this section we further comment on the reliability of our results at high redshifts 
and on the adopted extrapolations.
We also discuss some of the simplifications present in our method.
Finally, we place our results in the context of previous studies and
perform a brief comparison with the predictions of cosmological simulations.



 §.§ Reliability of our results at high redshifts


The current observational estimates of the GSMF used in our study reach to
z∼ 9 and the rate of redshift evolution of the SFMR
was also constrained at these high redshifts <cit.>.
However, the shape and evolution of the gas-phase MZR is essentially
unconstrained above z∼3.5 with the current observations (see Sec. <ref>)
and our results at higher redshifts need to be taken with a grain of salt.

 
The observational estimates of 
the total cosmic SFRD <cit.>
can serve as constraints to our model,
allowing to validate the combined assumptions about the SFMR and GSMF.
SFRD(z) recovered from our model shows a reasonable agreement
with the observations up to z≳8 for all variations in which
the low mass end slope of the GSMF for star forming galaxies does not evolve with redshift
(α_fix=-1.45).


However, there seems to be overall trend for α_fix to steepen with redshift 
(see Sec. <ref> and Fig. <ref>).
When we allow α_fix to evolve according to the fit presented in Fig. <ref>,
the cosmic SFRD(z) starts to significantly deviate from the observational estimates around z≳4
(see Sec. <ref>).
The cause of this discrepancy is not clear.


One of the assumptions which may be responsible for this problem
is that the low-mass end slope of the SFMR (a=0.83) does not evolve with redshift. 
As discussed in Sec. <ref>, it is difficult to constrain the evolution of a
with the current observations.
If the value of a would increase with redshift,
the amount of star formation occurring in low mass galaxies could be reduced and
this would alleviate the tension with the cosmic SFRD. 
Such evolution has been reported by several authors <cit.>,
while some other studies observe the reverse trend <cit.>.
Another argument in favor of the steepening a comes from the studies 
comparing the evolution of the stellar mass buildup resulting from 
the SFMR and as is observed from the evolution of the GSMF <cit.>.
The observationally inferred SFMR is known to lead to an overabundance of the low mass galaxies
with respect to the GSMF even when the high merger rate 
(shifting some of the excess galaxies to higher masses) is assumed for those galaxies.
<cit.> show that a∼1 at low masses 
(especially at the high redshift end of the redshift range considered in their study z≲2 ) 
helps to improve the consistency (although does not entirely solve the problem).


The inconsistency found in our variations with α_fix=α_fix(z)
may also arise from the assumed extrapolations to low stellar masses.
We extend the scaling relations and GSMF down to M_∗=10^6, which is much lower
than the stellar mass that can be probed by any high-redshift survey.
In principle, there is no guarantee that the slope of the GSMF and SFMR does not change
at stellar masses below the range probed by the current observations or that such a turnover
appears at high redshifts.
In fact, recent studies by <cit.> and <cit.> investigating 
the galaxy luminosity function (LF) in the extremely faint regime
with lensed galaxies find support for the presence of a flattening/turnover in the LF
at z∼6 and absolute UV luminosities M_UV≳ -15 mag
[
According to the mass - luminosity relation at z=6 from <cit.>,
M_UV∼ -15 mag and M_UV∼ -14 mag would correspond to
logM_∗∼6.5 and logM_∗∼6 respectively
and hence is close to the lower limit of the stellar mass of galaxies
considered in our study.
].
The upper limit on the UV emissivity at z∼5-6 found by <cit.>
also favors the presence of a turnover in the LF at M_UV≳ -14 mag.
However, the flattening/turnover at these low luminosities is not likely to solve
the discrepancy observed in our model.


Finally, we note that systematic trends in the IMF (e.g. with metallicity/redshift), if present,
could have a strong and non-straightforward effect on our results (see Sec. <ref>).




 §.§ Comparison with the Illustris cosmological simulations

While the comparison with the predictions of cosmological simulations
cannot settle the correctness of our results, it is certainly interesting
to confront the two methods, especially in the regimes not covered by observations.
Here we briefly compare our results with the predictions of the large-scale Illustris-1
cosmological hydrodynamical simulations <cit.> presented in Fig 1. in <cit.>,
leaving a more detailed comparison with the simulations for future studies.
The comparison is shown in Fig. <ref>.
The right panel reproduces Fig. 1 from <cit.> with additional
curves corresponding to the results from our model.
We show only three (out of six) lines presented in Fig. 1 in <cit.>
to improve the readability of the plot.
Those are the lines corresponding to the total SFRD (no cut
in ; green lines), the SFRD occurring below Z=0.0038 
(0.3 cut in  - =0.0127 in their paper; brown lines)
and the SFRD occurring below Z=0.0013 (0.1 cut in ; orange lines).
None of the curves corresponds to super-solar metallicity
(as  focus on long GRB host galaxies) 
and hence we cannot compare the distributions in that regime.
The comparison reveals that the closest match between the Illustris-1 results
and our model can be found for the low metallicity extreme variation
(the ratio stays within a factor of 2.5 up to z∼8).
The ratio changes only slightly if a different SFMR flattening is used in combination
with the PP04 MZR 
(e.g. no flattening in SFMR instead of the sharp flattening used in the low-Z extreme
variation; compare solid and dashed lines in the middle panel on the left 
in Fig. <ref>).
The best agreement is reached for the Z=0.0038 metallicity cut.
The SFRD occurring below Z=0.0013 in this variation is underpredicted
with respect to the Illustris-1 simulations at  z≲2 and overpredicted
at  z≳3, which suggests that the chemical evolution of the simulated
universe is slower than in our case.
In any case, at z∼8 the SFRD from our model shows a steep decrease - the
consequence of the evolution of the GSMF - not reflected in the results from simulations.
The Illustris-1 simulations predict higher star formation rate density
at z≲2 than our model and also more low metallicity star formation
up to z∼5 and above z∼8 for the moderate and high-Z extreme variations.


We note that the results of the comparison would likely change for the
updated IllustrisTNG simulations <cit.>, 
as one of the primary differences with respect to Illustris-1
is connected to the mass-metallicity relation from the simulations <cit.>.
<cit.> report the improved agreement with the observed MZR,
e.g. the flattening of the relation at high masses, which was not
reproduced in the previous version <cit.> is now accounted for.
Also, the slope of the simulated MZR has changed.
The normalization of the MZR found in IllustrisTNG continuously decreases
with redshift, with the average rate at z>2 of ≈ -0.064 dex per Δ z=1.
  We note that this rate is much lower than assumed in our model based on the 
   observed evolution between z=2.2 and z=3.5
  (see Sec. <ref> and <ref>).





 §.§ Our results in the context of previous studies
 



Several other studies used empirical scaling relations for star forming galaxies
to obtain the information about the distribution of the star formation
(or the formation of transients related to massive stars) across different metallicities.
Below we highlight the similarities and differences with respect to our method.
The details of the different methods are summarized in Tab. <ref> in appendix <ref>.
We note that none of those studies set finding the SFRD(Z,z) distribution
as its main goal and only a very limited comparison of our results with
the published material is possible. 

 <cit.> and <cit.> limited their analysis to low redshifts
(z≲0.3) and focused on application to long GRBs. 
Both studies take into account the intrinsic scatter in the SFMR and MZR.
<cit.> does not take into account the correlation between the SFR and metallicity
due to FMR, while <cit.> shows a variation
in which they use the FMR from <cit.> instead of combining the SFMR and MZR.
<cit.> consider two variations of the high-mass end slope of the SFMR -
relation with a downturn (also used by , see Fig. 1 therein)
and relation with flattening (similar to our sharp flattening case).
We note that such a downturn has not been observed in any of the subsequent studies
of the SFMR. This assumption may underestimate the star formation
happening at high metallicities.

 <cit.> and <cit.> use a slightly different approach than presented in this study:
they first obtain the metallicity distribution at each redshift (convolving MZR and GSMF)
and scale it with SFRD(z) to obtain the SFRD(Z,z)
 [ 
See also the study by <cit.> released during the revision process of this manuscript. 
We do not discuss their method here, but include it in the summary shown in Tab. <ref>.
].
Those authors constructed the SFRD(Z,z) that can be extended to 
high redshifts, albeit with various simplifications.
One of them is that the adopted mass metallicity relation assumes no scatter so that there is a unique
correspondence between mass and metallicity, and the redshift evolution only affects the normalization of the relation. 
Furthermore, both methods implicitly assume a constant specific star formation rate 
(i.e. SFR(M_∗)∝ M_∗; this is similar to our 'no flattening' case but with a different slope)
at all redshifts. 
They do not account for the FMR. 

<cit.> use MZR with the functional from <cit.> and redshift-dependent normalization,
adjusted to produce average metallicity of galaxies at z=0 of Z=0.03
or 0.016 in their high-end and low-end model respectively
(see Fig. 2 therein). 
<cit.> use a non-evolving GSMF and 
assume that the MZR is a single power law (without the observed flattening at high M_*)
and its normalization decreases with redshift (at z≳2
the decrease is much slower than in any variation considered in our study).

Note that the assumed lack of flattening in MZR (and possibly SFMR) does likely not affect the results of studies
that use only the low-metallicity part of the distribution (as in some application to long GRB e.g. ),
but may be important for other applications 
(e.g. the model by  has been recently used by to calculate 
the merger rates of double compact objects in , 
this was also the application of the model by ).

 Another example of using scaling relations to infer the information about metallicity 
comes from <cit.>. 
However, the authors give only the mass-averaged mean gas-phase metallicity
and not the entire distribution.
They use the mass-metallicity relation from <cit.> and GSMF from several studies
to estimate the mean at different redshifts, assuming that the relation given by <cit.>
can be extrapolated to z∼7.
The resulting mass-averaged mean metallicity is overplotted in Fig. <ref> 
showing our moderate variation. 
This mean is very close to the peak metallicity in the high metallicity extreme
(note that the MZR assumed in <cit.> is based on KK04 calibration,
as in the high-Z extreme case)
from our study at z≲3 and is significantly higher at higher redshifts.

 Using the results from <cit.>,
we construct the cumulative distribution of total star formation
below a certain metallicity for their method and plot it in Fig. <ref>
together with the ones from <cit.> (Fig. 2 therein, thick solid line[
The cumulative distribution for the case with SFMR that flattens at high masses 
discussed in <cit.> is not shown and hence we only compare with their
default case involving SFMR with a downturn at high masses.
]), <cit.> (Fig. 4 and 5 therein, solid lines) and resulting from our model.
It can be seen that all the curves fall within the broad range defined by the extreme cases
from our study (gray area).
The hatched region shows how much the result would change only due to changes in the SFMR
on the example of the high-Z extreme model variation.
The width of this range is much smaller than the width of the gray area,
which shows that it is mostly set by the differences in the MZR.
The results from <cit.> were obtained for the same metallicity calibration 
as assumed in the high-Z extreme case.
They reveal a noticeably greater offset towards low metallicities
from the high-Z extreme curve than indicated by the hatched region, 
especially if the correlation between the SFR and Z_O/H
is accounted for (as in our model).
This difference may be (partially) explained by the high-mass turnover
in the SFMR used by <cit.>, although as discussed earlier in this section,
there are also other differences in the assumptions
between our study and <cit.>.


Based on the data shown in the bottom panel of Fig. 2 from <cit.>, we conclude that
their low-end model would fall well within the gray area shown in Fig. <ref>.
However, the high-end model from <cit.> predicts hardly any z∼0
star formation below Z_O/H∼8.7 (see the middle panel in Fig. 2 therein) and therefore the 
cumulative distribution of the total star formation below certain Z_O/H
for that model would be shifted to noticeably higher metallicities with respect to our high-Z 
extreme shown in Fig. <ref>.



§ CONCLUSIONS

We combine empirical scaling relations for star forming galaxies
over wide range of masses and redshifts and other observational properties
of those galaxies to address the question
of (i) how was the cosmic star formation rate density distributed over metallicities
throughout the history of the Universe (SFRD(Z,z)) and
(ii) how uncertain is that distribution
given the currently unresolved problems
concerning the observationally inferred scaling relations.


To evaluate the uncertainty of the SFRD(Z,z) obtained with our method,
we use different versions of the scaling relations, covering the range of
possibilities present in the literature.
We use four different versions of the mass - metallicity relations
(based on different metallicity calibrations; 
with different normalizations, shapes and evolution with redshift),
vary the assumptions about the high mass end of the star formation - mass relation
(assuming either a single power law, broken power law or sharp flattening at the high mass end)
and about the low mass end of the galaxy stellar mass function
(assuming a fixed low mass end slope or a slope that steepens with redshift).
We take into account the scatter present in the relations,
the observed correlation between the SFR and metallicity 
(the fundamental metallicity relation)
and the distribution of gas metallicity within galaxies.
We extend the analysis up to z∼10, but we note that the MZR is not constrained by
observations above z∼3.5 and hence our results at higher redshifts rely on extrapolations.
To quantify the differences between the variations of our model, we compare the 
fractions of stellar mass formed at low (below 10% solar) and high (above solar) metallicity
at z≤0.5, 3 and 10. 
Our main results can be summarized as follows:

 
  * We identify two variations - the low and high metallicity extremes - that
maximize the fraction of stellar mass formed at low and high metallicity respectively.
Those extremes are reached by combining the PP04 MZR 
(based on the  O3N2 metallicity calibration)
with the SFMR with sharp flattening at high masses in the low-Z extreme case
and KK04 MZR 
(based on the  metallicity calibration)
with the SFMR with no flattening at high masses in the high-Z extreme.
The fraction of mass formed in stars at low (high) metallicity since z=3 
differs by ∼18% (∼26%) between the two extremes.
Those two variations also lead to different total stellar mass formed at each redshift
(with the higher stellar mass formed in the high-Z extreme).
This reduces the difference between the extremes if the stellar mass (instead of the fraction) 
formed at low metallicity is considered and
increases the difference between the amount of stellar mass formed at high metallicity.

  * We find that the variations with the low mass end slope of the GSMF that steepens with redshift
significantly overpredict the total star formation rate density at z≳4
with respect to other observational estimates. The cause of this discrepancy is not clear.
Those variations were not considered at z≳4.

  * The differences between the MZRs obtained for different metallicity calibrations
(in particular the differences in normalization) 
are the main cause of uncertainty in our results at z≲4.

  * We compare the results from our model with the local specific core-collapse supernovae (CCSN) rate
as a function of metallicity from <cit.>.
The uncertainties in the CCSN-location metallicity estimates, as well as in the
translation of the total SFR to the number of CCSN do not allow us to draw any 
strong conclusions from that comparison.
However, it might suggest that the variations with the sharp flattening 
at the high mass end of the SFMR underestimate the local high metallicity star formation.
The other variations of our model show a reasonably good agreement with the data.

Our model is publicly available and can be used in studies
focusing on the properties of populations of systems composed of stars and their remnants,
stellar-evolution related transients and their likely environments.
In particular, it can be applied to calculate the rates of those transients and to evaluate
their uncertainty due to the assumed distribution of the cosmic star formation rate density 
across metallicities and redshifts.
It can also serve for comparison for the cosmological simulations.





§ ACKNOWLEDGEMENTS

We would like to thank Leindert Boogaard, Søren Larsen, Svea Hernandez, Tereza Jeřábková, Jorryt Matthee, David Carton
and especially Jarle Brinchmann for valuable discussions. 
We also thank Or Graur, Zhiqiang Yan and the anonymous referee for their comments and Or Graur and Lucas Bignone
for sharing the data used in Fig. 11 and Fig. 12 respectively.
MC and GN acknowledge support from the Netherlands Organisation for Scientific Research (NWO).




mnras
 









§ METALLICITY - NOMENCLATURE AND SOLAR COMPOSITION

The term `metallicity` has been used variously in the literature,
depending on the object of interest and the method used.
It is generally understood as a measure of
the abundance of metals (elements heavier than helium) in a certain object or system
and defined as a fraction of total baryonic mass contained in metals:

    Z≡M_metals/M_baryons

although this definition is mostly applicable in theoretical studies.
Observationally metallicity is often expressed in terms of the of the oxygen to hydrogen 
abundance ratio:

    Z_O/H≡ 12 +  log (O/H ) ≡ 12 + log( n_O / n_H)

where n_i stands for the number density of either oxygen or hydrogen.
This definition is commonly used e.g. in interstellar medium studies, since
oxygen is the most abundant heavy element and it is relatively easy to observe
due to its strong optical lines. The common assumption is that the other elements
scale linearly with the measured Z_O/H maintaining the solar abundance ratios
(so that the conversion between the two measures has the form:
log(Z/Z_⊙) = Z_O/H - Z_O/H ⊙, where Z_⊙ and Z_O/H ⊙
stand for the solar metallicity in both measures).


However, the exact value of solar metallicity and solar composition is presently unknown. 
For instance, a widely used result obtained by <cit.>
based on 3D modeling of the solar atmosphere gives Z_⊙=0.0134 
and Z_O/H ⊙=8.69.
 Such low values of metallicity made it difficult
 to recover the helioseismically inferred information about the solar interior,
 such as sound speed profile, the location of the convective zone boundary and surface helium abundance
 <cit.>.
 In fact, the metallicity estimates based on helioseismology are systematically higher,
  with Z_O/H ⊙∼8.85 <cit.>.
 Using the abundances found by <cit.>, the resulting Z_⊙≈0.019.
 Similarly high values were obtained by the earlier solar atmosphere studies, 
 e.g. <cit.> found Z_O/H ⊙=8.83
and <cit.> found Z_O/H ⊙=8.93.
Asplund et al. (2009) recalculated Z_⊙ from those studied using the updated He abundances
 from helioseismology and found Z_⊙=0.017 and Z_⊙=0.0201 respectively.



§ INITIAL MASS FUNCTION - CORRECTIONS

The choice of a particular IMF impacts the estimates of both stellar masses and star formation rates
inferred from observations and as a result affects all of the galaxy scaling relations used in our model.
The results obtained assuming a certain IMF can be
translated to another IMF by shifting the mass and SFR estimates by a 
certain value (the change is in the same direction, but generally with different magnitude).
To correct for differences between the IMFs used in various studies 
(either <cit.>, <cit.> or <cit.>) 
we use the correction factors listed in table <ref>, unless explicitly stated otherwise.
These include the mass conversions as given in <cit.> (see eq. 2 therein)
and simple SFR correction calculated as the ratio of the total mass in both IMFs 
assuming the same number of massive stars 
<cit.>. 

 



§ METALLICITY OF STARS - DIFFERENCES BETWEEN THE VARIATIONS



 §.§ The metallicity at which the most of the stellar mass was formed


As discussed in Sec. <ref>, 
in almost all cases discussed in our study most of the stellar mass forms outside the assumed
low and high metallicity regimes (i.e. between 0.1  and ).
Here we look at the exact location of the peak of the distribution of mass formed in stars across metallicities
for different redshift ranges and discuss its dependence on the assumptions about the MZR, SFMR and GSMF.


Figure <ref> shows an example of such distribution.
We note that for all the variations considered in our study
such distributions can be described by a combination of a left-skewed normal 
distribution and a power-law tail at low Z_O/H.
When higher redshifts are included, the peak of the skewed-normal part shifts towards lower metallicities
 and the distribution broadens on the low metallicity side while the tail builds up.

 


The peak metallicity for different model variations is shown in Fig. <ref>.
It shows relatively little dependence on the assumed shape of the high mass part of the SFMR.
This can be seen by comparing at the location of the symbols with the same colour and different
shapes in Fig. <ref>.
The difference is ∼0.1 dex. 
Comparing points plotted with the same symbol and colour but
with and without edge in Fig. <ref> we conclude that
the effect of the low mass end slope of the GSMF is negligible.
To investigate the effect of the assumed mass metallicity relation 
(MZRs obtained for different metallicity calibrations),
we compare the locations of points plotted with the same symbol but with different colours.
It can be seen that for the T04, KK04 and M09 MZRs the difference is within 0.1 dex.
Those three relations have similar (high) normalizations at z∼0.
Additionally, both the M09 and T04 relations are relatively steep and evolve with redshift at a similar rate
and hence the green and blue points corresponding to those relations can be found close to each other or 
they even overlap in Fig. <ref>. The KK04 MZR is flatter than those two relations
(so that according to this MZR low mass galaxies have higher metallicities compared to M09 or T04)
and shows a milder redshift evolution and its peak metallicity is the highest.
The peak metallicity obtained for the PP04 MZR (with the low initial normalization) falls ≳0.25 dex 
below the other three MZRs.


 



 §.§ Low metallicity mass fraction
 
The fraction of mass formed in stars at low metallicities since z=0.5, 3 and 10 for different variations 
of the model is shown in Fig. <ref>.
The choice of the MZR (or metallicity calibration) has a decisive role in setting
the fraction of the low metallicity star formation 
(compare the same symbols plotted with different colours).
As expected, this fraction is the highest for PP04 (brown) and the lowest for
the KK04 MZR (red) at all redshifts.

The second important factor is the flattening at the high mass end of the SFMR.
For a given MZR (fixed colour in Fig. <ref>) the sharp flattening case (circles)
leads to the highest fraction of mass formed at low metallicity
 and this fraction is the lowest if the SFMR is a single power-law (`no flattening` case; squares).


The evolving low mass end slope of GSMF leads to a slight increase in the amount of mass
formed at low Z_O/H since z=3 with respect to the fixed slope case 
(as there are more low mass galaxies forming low metallicity stars).
This variation is not considered at higher redshift for the reasons discussed in sec. <ref>.

 



§ COMPARISON WITH CORE-COLLAPSE SUPERNOVAE RATES

 
 

 §.§ Cosmic SFRD vs CCSN rate density

It is commonly assumed that the CCSN rate is proportional to the SFR and the
constant of proportionality depends solely on the IMF and can be
estimated as a number of stars that explode as CCSNe per unit mass 
<cit.>.
This calculation requires assuming a certain initial mass range for the progenitors of CCSNe.
However, the exact limits of this mass range are not known 
<cit.>.
Furthermore, those limits depend on the stellar metallicity <cit.>
and hence likely change throughout the cosmic history.
The expected CCSN progenitor mass range is further modified by binary interactions and stellar mergers
<cit.>.
As argued by <cit.>, over 70% of the most massive stars are born in interacting binary systems
and ∼20 - 30% will experience a merger with their companion.
Stellar interactions and mergers can act both to reduce the number of CCSN expected from a certain stellar population
(e.g. by mergers of two CCSNe progenitors) and to increase that number 
(e.g. by mergers of the more common low mass stars
that result in a CCSN progenitor). 
<cit.> find that the overall effect of accounting for the binaries
is to increase the number of the CCSN expected to form from a stellar population of the same mass
(they find ∼14% increase, although the exact number is model dependent).
Another complication added by binaries is the fact that 
for the same amount of star formation, stellar population containing binaries would lead to a
higher UV luminosity than the one consisting only of single stars and therefore would lead to an overestimate
in the SFR measured using the common UV-light based tracers 
<cit.>.
In other words, the actual number of stars born as CCSN progenitors would be smaller than what results from 
the inferred SFR, which would lead to an overestimate of the number of expected CCSN.
The combined effect of the presence of binaries in the observed stellar population on the 
CCSN formation efficiency is not clear.
The comparison is further complicated by the fact that the
observational CCSN rate estimates are not always consistent with each other, 
which may be partially attributed to different methods used to correct for the number of CCSN
missed due to obscuration by dust <cit.>.
 


 §.§ Central metallicity vs metallicity at the supernova location


<cit.> suggest that the true metallicity at the CCSN location should be contained within
the range reaching metallicity 0.1 dex lower than the central value.
However, the <cit.> metallicity calibration used in <cit.> probes a relatively small metallicity range 
and as shown in <cit.>, leads to relatively flat metallicity gradient estimates.
For instance, <cit.> find the typical metallicity gradient of around -0.07 dex/r_e 
using the <cit.> calibration, -0.11 dex/r_e with the PP04 calibration and 
-0.14 dex/r_e using the theoretical calibration from <cit.>.
Hence, based on the metallicity gradient studies and taking this into account
the fact that CCSN are found to explode within 2 r_e from the centre,
<cit.>, 
one can expect to find also larger offsets between the central and local metallicity than 0.1 dex proposed by 
<cit.> (especially when using one of the theoretical metallicity calibrations).
<cit.> find a median offset between the nuclear and explosion site metallicity for their stripped envelope
SN sample of ∼0.08 dex in the PP04 scale, although with a substantial scatter (see Fig. 12 therein). 
However, the galaxies considered in this study contain objects at higher redshift and with smaller angular sizes
than those used in <cit.>. Hence, the nuclear metallicity in those galaxies was measured
with a larger covering fraction (averaging information from a larger portion of the galaxy) 
and can be expected to be closer to the metallicity at the SN explosion site.
<cit.> compared the local metallicities of part of their sample of stripped envelope supernovae with the
central metallicities obtained from the nuclear SDSS spectra using the PP04 metallicity scale and found the average 
difference of 0.13 dex and in the extreme cases reaching up to 0.24 dex.
 


 §.§ Specific CCSN rate vs metallicity - dependence on the scatter in metallicity


The shape of the dependence of the specific CCSN rate on metallicity is a consequence
of the scaling relations in star forming galaxies.
The specific star formation rate (SFR/M_∗) decreases with mass, so that there are fewer CCSN per stellar mass
in more massive host galaxies. The mass - metallicity relation tells that with the increasing stellar mass 
(and hence decreasing specific SFR or specific CCSN rate) metallicity increases up to a certain Z_O/H value,
at which the MZR flattens, hence the slope of the specific CCSN rate - metallicity relation steepens around the 
metallicity at which the MZR begins to flatten. 
The presence of the scatter around the MZR pushes a fraction of the star formation
into higher and lower metallicity than what results from the MZR, 
but the anticorrelation between the SFR and metallicity 
(resulting from the fundamental metallicity relation)
reduces that `smearing` effect of the scatter and acts to further steepen the high-metallicity slope.
The gray solid line in Fig. <ref> shows the specific CCSN rate - metallicity
relation when only the scatter in the MZR (σ_0) is taken into account.
It can be seen that the high-metallicity slope of the gray curve is steeper than in the case of the blue solid curve,
which corresponds to the same model variation, but taking into account the fact that the
stars within galaxies can form with a range of metallicities (σ_∇ O/H=0.14 dex).
Hence, the assumed width of this distribution has a noticeable effect on the high metallicity end slope
and on the location of the turnover in the specific CCSN rate - metallicity dependence.

For instance, the power law slope of the high metallicity Z_O/H≳ 8.5 part of this relation
is ∼ -1 for σ_∇ O/H=0.14 dex, ∼ -1.4 for σ_∇ O/H=0.1 dex
and for ∼ -0.65 for σ_∇ O/H=0.2 dex.  
 
 


§ CLOSER LOOK AT THE SECONDARY ASSUMPTIONS
 


 §.§ The distribution of the star formation over metallicities within galaxies

As discussed in Sec. <ref>, there is a range of metallicities at which
stars form within their host galaxies. 


Our treatment of the distribution of metallicity of the star formation within galaxies
is very simplified. We assume that it can be described as a normal distribution centered
at the metallicity resulting from the MZR 
(taking into account the anticorrelation between the metallicity and the SFR) and with a
dispersion σ_∇ O/H=0.14 dex for all galaxies and redshifts and 
irrespective of the choice of the MZR (metallicity calibration).

However, the range of metallicities found in the HII regions (setting the width of the distribution),
as well as the value of the characteristic metallicity gradient depend
on the assumed metallicity calibration <cit.>.
The value of σ_∇ O/H that we use is based on the averaged results obtained using
the <cit.> and the <cit.> metallicity calibrations and could be higher if different 
(theoretical) metallicity calibrations were used.


On the other hand, the offsets between
the metallicities measured at the location of the CCSN and the central metallicity of
its host galaxies (which also indicate the range of metallicities at which the stars form
within galaxies) reported in the studies using the empirical metallicity calibrations mentioned above
are smaller than 2×σ_∇ O/H (see appendix <ref> and references therein).


Furthermore, not all galaxies reveal a single slope negative metallicity gradient.
<cit.> studied galaxies at intermediate redshifts and found a significant fraction
of galaxies showing gradients that are consistent with being flat and
several cases of positive gradients, which might suggest that the simple common negative metallicity gradient
found in the galaxies in local Universe may not be relevant at higher redshifts.


Finally, the detailed studies of metallicity gradients 
within galaxies (on which we based our assumptions) are limited to the relatively local Universe.
It is thus relevant to ask how much our results depend on this particular assumption.


  §.§.§ The importance for our results

To asses the importance of the assumptions about the distribution of the star formation over metallicities
within galaxies on our results, we vary the value of σ_∇ O/H in the range 0.1-0.2 dex.
The effect on the results (summarized in Table <ref>) is minor:
the difference in the total stellar mass formed and the fraction of the low metallicity star formation
is negligible and the fraction of the high metallicity star formation differs by less than 3%
between the calculations with σ_∇ O/H=0.1 and 0.2 dex for all model variations.
The distribution of mass formed in stars at different metallicities since a certain redshift 
(see Fig. <ref>) slightly broadens as σ_∇ O/H increases,
however this broadening only affects the low metallicity (Z<0.1) tail of the distribution 
at high redshifts z>5.
At the same time the peak metallicity shifts to lower values 
(by less than 0.03 dex for the distribution at z<0.5 and less than 0.09 dex at z<10 for all model variations), 
which reduces the impact of the broadening on the fraction of the high metallicity star formation.
However, as discussed in Sec. <ref>, changing σ_∇ O/H noticeably affects
the slope of the high metallicity end of the specific CCSN rate - metallicity relation
(see Fig. <ref> and Sec. <ref>).



 §.§ Contribution from passive galaxies


In this study we use the scaling relations and galaxy mass functions obtained for star forming galaxies.
However, various criteria to select those galaxies are used in the literature and
a fraction of galaxies classified as passive may still have some residual star formation.
The fraction of passive galaxies increases towards lower redshifts <cit.>
and hence their contribution to the global star formation budget can be expected to be
the highest in the local Universe.
<cit.> use the SDSS DR7 release data to study the star formation - mass relation of
the local galaxies without imposing any cuts on the star formation within galaxies in their sample, 
i.e. including quenched galaxies.
They find that a low level of star formation is still present within the passive population 
(although they note that in many cases the measured SFR may be an upper limit)
and show that while the majority of mass in the local Universe reside in those galaxies, 
they do not contribute significantly to the ongoing star formation.


  §.§.§ Importance for our results

To obtain a rough estimate of the upper limit on the local star formation rate density
contributed by passive galaxies,
we construct the `upper limit star formation - mass relation` for passive galaxies by 
fitting a line connecting the two low SFR peaks in the SFR - stellar mass plane shown in Fig. 4 in <cit.>.
Combining it with the local galaxy stellar mass function for the passive (red)
population of galaxies from <cit.> and integrating over the stellar masses between
10^8 - 10^12, we find that the contribution of passive galaxies to the local SFRD
is more than a factor of 80 smaller than the SFRD obtained for the star forming galaxies in that mass range
in the moderate variation of our model.
Most of this residual star formation occurs in massive galaxies and hence adds to the high metallicity
tail of the SFRD(Z,z) distribution.
Assuming that all the star formation from passive galaxies adds to the star formation at Z>
and their contribution to the global SFRD is at the level corresponding to the estimated upper limit,
the fraction of mass formed in stars at high metallicity at z<0.5 in the moderate variation of our model
would increase by ≲1%.
Taking into account the fact that the fraction of passive galaxies decreases with redshift, 
only some of them show some residual star formation and the above estimate used the upper limit on the
level of that residual star formation rate, we conclude that the omitted star formation in passive
galaxies has negligible effect on our results.



§ SFRD(Z,Z) FROM SCALING RELATIONS - OTHER STUDIES


Table <ref> summarizes the assumptions made by several different authors using empirical
(or mixed empirical and theoretical/simulated) 
scaling relations to infer SFRD(Z,z).

 





