The GREATS H$\beta$+[OIII] Luminosity Function and Galaxy Properties at $\mathbf{z\sim8}$: Walking the Way of JWST

TL;DR

This study derives the [OIII]+Hβ luminosity function at z~8 using combined Spitzer and HST data, providing insights into galaxy properties and predicting JWST spectroscopic yields, highlighting the need for deep pre-imaging.

Contribution

First to establish the [OIII]+Hβ luminosity function at z~8 using combined photometric data and photoionization models, informing future JWST observations.

Findings

[OIII]+Hβ luminosity function is higher at z~8 than at lower redshifts.

Strong correlation between [OIII]+Hβ and UV luminosity at z~8.

Current legacy fields are too shallow for optimal JWST spectroscopy at z~8.

Abstract

The James Webb Space Telescope will allow to spectroscopically study an unprecedented number of galaxies deep into the reionization era, notably by detecting [OIII] and H$\beta$ nebular emission lines. To efficiently prepare such observations, we photometrically select a large sample of galaxies at $z\sim8$ and study their rest-frame optical emission lines. Combining data from the GOODS Re-ionization Era wide-Area Treasury from Spitzer (GREATS) survey and from HST, we perform spectral energy distribution (SED) fitting, using synthetic SEDs from a large grid of photoionization models. The deep Spitzer/IRAC data combined with our models exploring a large parameter space enables to constrain the [OIII]+H$\beta$ fluxes and equivalent widths for our sample, as well as the average physical properties of $z\sim8$ galaxies, such as the ionizing photon production efficiency with $\log(\xi_\mathrm{ion}/\mathrm{erg}^{-1}\hspace{1mm}\mathrm{Hz})\geq25.77$. We find a relatively tight correlation between the [OIII]+H$\beta$ and UV luminosity, which we use to derive for the first time the [OIII]+H$\beta$ luminosity function (LF) at $z\sim8$. The $z\sim8$ [OIII]+H$\beta$ LF is higher at all luminosities compared to lower redshift, as opposed to the UV LF, due to an increase of the [OIII]+H$\beta$ luminosity at a given UV luminosity from $z\sim3$ to $z\sim8$. Finally, using the [OIII]+H$\beta$ LF, we make predictions for JWST/NIRSpec number counts of $z\sim8$ galaxies. We find that the current wide-area extragalactic legacy fields are too shallow to use JWST at maximal efficiency for $z\sim8$ spectroscopy even at 1hr depth and JWST pre-imaging to $\gtrsim30$ mag will be required.

Equations

$$\log (\mathrm{L}([\text{O}\text{textsc}\mathrm{iii}]\lambda \lambda 4959,5007+\mathrm{H}\beta )/\mathrm{erg}{\mathrm{s}}^{-1})=0.86\pm 0.12\times$$

$$\log ({\mathrm{L}}_{\mathrm{UV}}/\mathrm{erg}{\mathrm{s}}^{-1})+33.92_{-1.27}^{+1.23}$$

Full text

††pubyear: 2018††pagerange: The GREATS H$\beta$+[O III] Luminosity Function and Galaxy Properties at $\mathbf{z\sim 8}$: Walking the Way of JWST–The GREATS H$\beta$+[O III] Luminosity Function and Galaxy Properties at $\mathbf{z\sim 8}$: Walking the Way of JWST

1 Introduction

Great progress has been made over the last two decades in our study of early galaxy mass assembly and the evolution of the cosmic star-formation rate density at $z\geq 4$ (e.g., Madau & Dickinson, 2014; Duncan et al., 2014; Salmon et al., 2015). However, so far, these studies mostly rely on the analysis of broad-band photometry only, given that current facilities only provide access to the faint, rest-UV emission lines in a small number of bright galaxies (e.g., Stark et al., 2017). While a wealth of photometric data are now publicly available, photometric studies can suffer from several caveats. The selection of high-redshift galaxies relies on the Lyman break (i.e., dropout selection; Steidel et al., 1996) that can possibly miss a significant fraction of galaxies (e.g., Inami et al., 2017), and the derivation of most of the galaxy physical properties relies on spectral energy distribution (SED) fitting that is affected by several degeneracies and strongly depends on assumptions (e.g., star formation history, metallicity, dust extinction curve; Finlator et al., 2007; Yabe et al., 2009; De Barros et al., 2014). Furthermore, the photometry can be contaminated by strong nebular emission lines that are ubiquitous at high-redshift (e.g., Chary et al., 2005; Schaerer & De Barros, 2010; Shim et al., 2011; Stark et al., 2013; Labbé et al., 2013; Smit et al., 2014; Shivaei et al., 2015; Faisst et al., 2016; Mármol-Queraltó et al., 2016; Rasappu et al., 2016). To account for their impact, either empirical (e.g., Schaerer & De Barros, 2009) or dedicated photoionization modeling (e.g., Zackrisson et al., 2001) have been used. However, direct observational access to these emission lines through spectroscopy will have to await the advent of the James Webb Space Telescope (JWST).

At lower redshift, where most ultraviolet, optical, and near-infrared emission lines can be observed either from the ground or space, lines are efficiently used to determine instantaneous star-formation rates (SFR) and specific star-formation rates (sSFR=SFR/M*⋆*; Kauffmann et al., 2004), to derive gas-phase element abundances (e.g., Tremonti et al., 2004), to accurately derive the dust extinction thanks to the Balmer decrement (e.g., Domínguez et al., 2013; Reddy et al., 2015), and to determine the main source of ionizing photons (star formation or AGN) with the BPT diagram (Baldwin et al., 1981). While direct observations of optical and near-infrared lines are out of reach for $z>4$ galaxies until the launch of JWST, one can take advantage of the impact of nebular emission lines on the photometry to probe these lines indirectly and to uniquely reveal some ISM properties of high-redshift galaxies.

Emission lines can also be useful to derive very accurate photometric redshifts. Smit et al. (2015) exploit extremely blue Spitzer/IRAC colors to identify $6.6\leq z\leq 6.9$ galaxies, for which [O iii]$\lambda\lambda 4959,5007$+H$\beta$ lines are expected to fall in the 3.6µm band while the 4.5µm band is free of line contamination. A similar technique applied to $z>7.1$ galaxies has led to the reliable selection and subsequent spectroscopic confirmation of some of the most distant Lyman-$\alpha$ emitters to date (Roberts-Borsani et al., 2016; Oesch et al., 2015; Zitrin et al., 2015; Stark et al., 2017). SFR and sSFR can be derived from H$\alpha$ emission for galaxies at $z\sim 4$ where the H$\alpha$ line is found in the IRAC 3.6µm channel while the IRAC 4.5µm channel is free from line contamination (Shim et al., 2011; Stark et al., 2013; Mármol-Queraltó et al., 2016). Furthermore, the H$\alpha$ luminosities have also been used to derive the ionizing photon production efficiency ($\xi_{\mathrm{ion}}$, defined as the production rate of ionizing photons per unit luminosity in the UV-continuum; Bouwens et al., 2016) at $z\sim 4$. The comparison between uncorrected SFR derived from emission lines and SFR derived from UV+IR shed light on the relative stellar to nebular attenuation (Shivaei et al., 2015; De Barros et al., 2016).

In this paper, we adopt the same approach as these latter studies: we use the impact of nebular emission on the broad-band photometry to indirectly derive emission line fluxes and EWs to gain insight into high-redshift galaxy properties. Specifically, we use a sample of photometrically selected $z\sim 8$ galaxies, for which the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ lines fall in the IRAC 4.5µm channel while we do not expect strong lines in the IRAC 3.6µm channel. We take advantage of the Spitzer ultra-deep survey covering CANDELS/GOODS South and North fields, the GOODS Re-ionization Era wide-Area Treasury from Spitzer (GREATS, Labbé et al. 2019, in prep) survey, providing the best constraints on $z\sim 8$ IRAC colors to date. To derive the line fluxes as accurately as possible and account for most of the uncertainties, we use synthetic SEDs produced with dedicated photoionization modeling to fit the observed $z\sim 8$ SEDs. Our aim is to derive the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ luminosity function (LF) at $z\sim 8$ to prepare efficient JWST observations in the future.

The paper is structured as follows: The photometric data and selection procedure are described in Sec. 2. We provide a description of our photoionization grid in Sec. 3, and Sec. 4 gives the SED fitting method. In Sec. 5, we present the constraints that we obtain on physical properties of $z\sim 8$ galaxies. The resulting [O iii]$\lambda\lambda 4959,5007$+H$\beta$ luminosity function is shown and discussed in Sec. 6. We summarize our conclusions in Sec. 7.

We adopt a $\Lambda$-CDM cosmological model with $\mathrm{H}_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{\mathrm{m}}=0.3$ and $\Omega_{\Lambda}=0.7$. All magnitudes are expressed in the AB system (Oke & Gunn, 1983).

2 Data and Sample

The input sample used for this work is based on the Lyman break galaxy (LBG) catalogs from Bouwens et al. (2015). These are compiled from all the prime extragalactic legacy fields, including the Hubble Ultra Deep Field (HUDF; Ellis et al., 2013; Illingworth et al., 2013) and its parallel fields, as well as all five CANDELS fields (Grogin et al., 2011; Koekemoer et al., 2011).

In addition to deep $HST$ near-infrared imaging, all these fields have extensive $Spitzer$/IRAC coverage. We have reduced and combined all the IRAC 3.6 and 4.5 $\micron$ data available in each field. In particular, we include the complete data from the GREATS survey (Labbe et al. 2019, in prep). GREATS builds on the vast amount of archival data in the two GOODS fields (Giavalisco et al., 2004) and brings the IRAC 3.6µm and 4.5µm coverage to a near-homogeneous depth of 200-250 hr, corresponding to 5$\sigma$ sensitivities of 26.8-27.1 mag, over $\sim 200$ arcmin2. This is very well matched to the $H_{160}$-band detection limits, allowing us to detect the rest-frame optical light of nearly all the galaxies identified in the $HST$ data. For a complete description of the $HST$ and $Spitzer$ dataset we refer the reader to Bouwens et al. (2015) and Stefanon et al. (2019, in prep).

Given the much wider point-spread function (PSF) of the IRAC data compared to $HST$, special care is required to derive accurate photometry. We use a custom-made software tool mophongo, developed and updated over the last few years (e.g., Labbé et al., 2010, 2015). In short, starting from the $HST$ F160W image, mophongo uses position-dependent $HST$-to-IRAC PSF kernels to fit and subtract all the neighboring galaxies in a 21″ region around a source of interest, before measuring its flux density in a 2″ aperture.

After discarding objects with IRAC fluxes highly contaminated by neighbours ($<10\%$), we ran the EAZY code (Brammer et al., 2008) on the full photometric catalog, allowing the photometric redshift to vary from 0 to 10, and deriving the redshift probability distribution for each object in the GREATS data. We took care to exclude IRAC 3.6µm and IRAC 4.5µm bands from the SED fit, since emission lines at high redshift can affect the photometry in these bands (e.g., Smit et al., 2014) and we want to avoid to be biased toward objects with large EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$). We want to focus on $z\sim 8$ galaxies for which the IRAC 3.6µm$-$IRAC 4.5µm color provides constraints on the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ flux and these lines have completely entered the Spitzer/IRAC 4.5µm channel at $z\geq 7.11$ 111The [O iii]+H$\beta$ lines are out of IRAC 4.5µm at $z>9.05$ and after applying our selection, only one galaxy has a redshift above this limit.. Therefore, the selection criterion was defined as $p(z\geq 7.11)\geq 0.68$, with $p(z\geq 7.11)$ the probability for a galaxy to have a redshift $z\geq 7.11$. Additionally to this criterion, we select a subsample of galaxies with at least one detection with $S/N\geq 3$ in either IRAC 3.6µm or IRAC 4.5µm  allowing us to derive at least a strong upper or lower limit on the $(3.6-4.5)$µm color for each of these galaxies. In this subsample ($N=76$), 16 galaxies are detected at 4.5µm only and 10 at 3.6µm only, potentially slightly biasing this subsample toward larger EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$), with the detection at 4.5µm being consistent with an increase in flux due to [O iii]+H$\beta$ lines, while the stellar continuum could remain undetected at 3.6µm. We show in Fig. 1 the F125W magnitude distribution for the galaxy selected based on their photometric redshifts ($N=135$) and the same distribution after applying the Spitzer/IRAC detection criterion. The F125W band probes the UV rest-frame emission of $z\sim 8$ galaxies ($\sim 1500$Å). $\sim 50\%$ of the $z\sim 8$ galaxies are detected with $S/N\geq 3$ in at least one Spitzer/IRAC band, mostly the brightest with F125W$<27.5$. While we apply a threshold detection in IRAC to select our subsample, we use all available photometry including bands with low S/N (<3) as well as non detections to perform the SED fitting. We use the subsample ($N=76$) with $S/N(3.6\mu\mathrm{m}\lor 4.5\mu\mathrm{m})\geq 3$ to constrain the $z\sim 8$ galaxy physical properties (Sec. 5) and the entire photometric sample ($N=135$) to derive the [O iii]+H$\beta$ LF (Sec. 6).

3 Photoionization models

Several works have used photoionization models to study or predict nebular emission line properties of high-redshift galaxies (e.g., Zackrisson et al., 2011; Jaskot & Ravindranath, 2016; Steidel et al., 2016; Nakajima et al., 2018; Berg et al., 2018). Since we do not have access with the current facilities to the optical/NIR nebular emission lines for galaxies at $z\sim 8$, the ISM physical conditions at these high redshifts are largely unknown. Therefore we created a grid of photoionization models with a large parameter space to encompass the plausible stellar and ISM physical properties at $z\sim 8$

We used the latest release of the Cloudy photoionization code (C17, Ferland et al., 2017) to build our grid of models. We chose SEDs from the latest BPASSv2.1 models (Eldridge et al., 2017) as input, which account for stellar binaries and stellar rotation effects. Indeed, several recent observations point out that high-redshift galaxies can exhibit UV emission lines requiring hard ionizing photons (e.g., C iii]$\lambda\lambda 1907,1909$, C iv$\lambda\lambda 1548,1550$; Stark et al., 2014; Amorín et al., 2017; Vanzella et al., 2017), and there are mounting evidences that these UV lines are due to star formation since they are spatially associated with star forming regions (Smit et al., 2017). These kind of strong emission lines can be reproduced by including stellar rotation and/or stellar binaries in the stellar models (e.g., Eldridge et al., 2008) or using stellar templates updated with recent UV spectral libraries and stellar evolutionary tracks as in the latest Charlot & Bruzual single star models (e.g., Gutkin et al., 2016). Regarding the ISM properties, we adopt the same interstellar abundances and depletion factors of metals on to dust grains, and dust properties as Gutkin et al. (2016). These authors show that these modeling assumptions span a range that can reproduce most of the observed UV and optical emission lines at low- and high-redshift (Stark et al., 2014, 2015a, 2015b, 2017; Chevallard & Charlot, 2016). While Gutkin et al. (2016) use different stellar population synthesis model than used here, namely an updated version of the Bruzual & Charlot (2003) stellar population synthesis model, a comparison of these two SPS models show that they provide similar results in interpreting stellar and nebular emissions of local massive star-clusters (Wofford et al., 2016). For the grid used in this work, we use stellar metallicities from $Z=0.001$ to $Z=0.008$ with an initial mass function (IMF) index of $-2.35$ and an upper mass cutoff of 300M*⊙. For each stellar metallicity, for simplicity, we assume the same gas-phase metallicity. We also explore a range of C/O abundance ratio (from $\log\mathrm{C}/\mathrm{O}=-1.0$ to -0.4, consistent with the observations of Amorín et al., 2017), three different values of dust-to-metal ratios ($\xi_{d}=0.1,0.3,0.5$; Gutkin et al., 2016), and a range of hydrogen gas densities ($10^{2}$ to $10^{3}$ cm-3*). We assume no leakage of ionizing photons. For each set of parameters and each stellar age, we built SEDs, adding to the pure stellar SEDs from BPASS the nebular emission lines and nebular continuum as computed in Cloudy.

4 SED fitting

We use the spectral energy distributions created with Cloudy to perform SED fitting of each individual galaxy in our $z\sim 8$ sample. We use a modified version of the SED fitting code Hyperz (Bolzonella et al., 2000) allowing us to apply two different attenuation curves to the stellar and nebular components of the SED. We apply a Calzetti attenuation curve (Calzetti et al., 2000) to the stellar component and a Cardelli attenuation curve (Cardelli et al., 1989) to the nebular component, adopting $E(B-V)_{\mathrm{gas}}=E(B-V)_{\mathrm{\star}}$ for simplicity. Studies of galaxy samples at $z\sim 2$ show that the ratio $E(B-V)_{\mathrm{gas}}/E(B-V)_{\mathrm{\star}}$ is affected by a large scatter and is increasing with increasing $E(B-V)_{\mathrm{\star}}$ and SFR (Reddy et al., 2015; Theios et al., 2019). Since $z\sim 8$ galaxies exhibit blue UV $\beta$ slopes indicating low dust extinction (e.g., Bouwens et al., 2014), we do not expect dust to have a large impact on our results. Nevertheless we allow $E(B-V)$ to vary from 0.0 to 0.2 in our SED fitting procedure. We assume a constant star formation history. While the choice of the SFH has an impact on the derived physical parameters, this effect is alleviated for young ages ($<100$Myr; De Barros et al., 2014). Since the oldest age allowed by the cosmological model adopted in this work is 730 Myr and due to the relatively young best-fit ages found for our sample (Sec. 5), the assumed SFH has a limited impact on our final results.

Minimization of $\chi^{2}$ over the entire parameter space yields the best-fit SED. For each physical parameter of interest, we derive the median of the marginalized likelihood, and its associated uncertainties. Derived physical parameters include age of the stellar population, stellar mass, SFR, sSFR, as well as observed (i.e., attenuated by dust) [O iii]$\lambda\lambda 4959,5007$ and H$\beta$ emission line fluxes and EWs, and ISM physical properties (e.g., ionization parameter). All physical properties used in this work such as EW([O iii]+H$\beta$) and L([O iii]+H$\beta$) are SED derived, and correspond to the median of the marginalized likelihood, except stated otherwise.

The median uncertainty regarding the EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$) for the entire sample is $+0.34/-0.37$ dex. We also split the sample in four bins of F125W magnitude to emphasize the reliability of the constraints depending on the UV luminosity. For $\mathrm{F125W}<26$, $26<\mathrm{F125W}<26.5$, $26.5<\mathrm{F125W}<27$, and $\mathrm{F125W}>27$, the median uncertainties are $+0.27/-0.32$ dex, $+0.30/-0.32$ dex, $+0.38/-0.45$ dex, and $+0.39/-0.52$ dex, respectively. We describe how we account for those relatively large uncertainties in the [O iii]+H$\beta$ luminosity function (LF) derivation in Sec. 6.1.

5 Constraints on the ISM and Physical Properties

5.1 Predictions for the (3.6-4.5)µm Color

We show in the top left panel of Fig. 2 the range of (3.6-4.5)µm color which is spanned by our grid of models for $Z=0.004$. At an age of 1 Myr, (3.6-4.5)µm can vary by 2 magnitudes, nebular emission (mainly [O iii]$\lambda\lambda 4959,5007$+H$\beta$ lines) boosting the flux in the IRAC 4.5µm channel and producing IRAC 3.6µm$-$IRAC 4.5µm color redder by $\sim 1.5$ magnitude in comparison with the color expected from a pure stellar template. However, nebular emission can also have the opposite effect for low ionization parameters, producing a bluer color than expected for pure stellar emission up to $\sim 0.5$ magnitude. This effect is due to the relation between the ionization parameter and the [O iii]/[O ii]$\lambda 3727$ ratio: $\log([\textrm{O}~{}\textsc{iii}]/[\textrm{O}~{}\textsc{ii}]\lambda 3727)$ is increasing with higher ionization parameter (depending on the metallicity, Kewley & Dopita, 2002), becoming larger than 1 at $Z=0.05Z_{\odot}$ for $\log\mathrm{U}\gtrsim-3.0$. However, the number of galaxies for which the IRAC 3.6µm$-$IRAC 4.5µm color is best fitted with a very low ionization parameter ($\log\mathrm{U}=-4.0$) is small (25%). Furthermore, photometric redshift uncertainties can also account for $z>7.1$ blue colors with the Balmer jump starting to enter the IRAC 3.6µm channel at $z>8$.

We show the EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$) evolution with age in Fig. 2 (bottom left panel) and examples of SEDs for our final sample for a range of F125W magnitude in Fig. 3.

5.2 Stellar Metallicity and ISM Physical Properties

The models required to reproduce the SEDs, mainly the IRAC colors, give insight in the ISM physical properties at $z\sim 8$. Some parameters, like the C/O ratio or the hydrogen gas density, have little to no impact on the EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$) which is the spectral feature with the largest effect on the IRAC colors, and therefore providing the main constraints on the ISM physical conditions. In our grid of models, the (3.6-4.5)µm is mostly defined by the stellar metallicity, the ionization parameter, and the age of the stellar population. For our sample, we find a median stellar metallicity of $Z_{\star}=0.004_{-0.002}^{+0.004}$, a median ionization parameter $\log\mathrm{U}=-3.0\pm 1.0$, and a median $\log(\mathrm{age}/\mathrm{yr})=7.2^{+0.9}_{-0.6}$. The constraints on dust extinction are mostly coming from the fit of the UV $\beta$ slope and we find a median extinction $A_{V}=0.4\pm 0.2$.

As noted previously, we do not assume any ISM properties of $z\sim 8$ galaxies but compare their photometry with a grid of photoionization models that we consider to encompass the plausible $z\sim 8$ properties. The constraints on the stellar metallicity and the ionization parameter are driven by the number of ionizing photons available to interact with the gas and the gas-phase Oxygen abundance. The ionizing photon flux is increasing with decreasing metallicity (e.g., Stanway et al., 2016) and with increasing ionization parameter, but the gas-phase Oxygen abundance is decreasing with decreasing metallicity, and also decreasing with increasing dust-to-metal ratio (more Oxygen depletes into the ISM dust-phase; Gutkin et al., 2016). Therefore the ISM parameter derivation suffers from several degeneracies: to reproduce (3.6-4.5)µm colors produced by strong [O iii]$\lambda\lambda 4959,5007$+H$\beta$ emission lines ($(3.6-4.5)\micron>0$, Fig. 2) a large range of metallicities and ionization parameters is allowed, as long as there is the right balance between the number of ionizing photons and the Oxygen abundance. However, the parameter space allowing this balance is smaller for low metallicities due to large ionizing photon production but low Oxygen abundance. The same is true for high metallicities due to lower ionizing photon production and high Oxygen abundance. Then the median metallicity found for our sample only reflects that for $Z=0.004$ there is a larger parameter space in terms of ionization parameter, age, and dust-to-metal ratio allowing to reproduce the observed (3.6-4.5)µm colors. Indeed, we found that $\sim 40\%$ of our sample has a best-fit SED with $Z\leq 0.002$, while the median of the marginalized likelihood for the metallicity is $Z\geq 0.004$ for the entire sample.

The choice of an IMF upper mass cutoff at 300M*⊙* has a negligible impact on our results since a cutoff of 100M*⊙* changes the [O iii]$\lambda\lambda 4959,5007$ flux by $\sim 10\%$ (1-2% for H$\beta$ and [O ii]) which is small compared to the typical uncertainties affecting EWs and line luminosities (Sec. 4). Changing the assumed dust attenuation curve from a Calzetti to an SMC curve (Prevot et al., 1984; Bouchet et al., 1985) has also little to no impact on the overall derived properties, except for dust attenuation.

The $z\sim 8$ subsample with $S/N(3.6\mu\mathrm{m}\lor 4.5\mu\mathrm{m})\geq 3$ has a median stellar mass $\log(\mathrm{M}_{\star}/\mathrm{M}_{\odot})=8.62^{+0.43}_{-0.39}$ and a median SFR $\log(\mathrm{SFR}/\mathrm{M}_{\odot}{\rm yr}^{-1})=1.26^{+0.42}_{-0.30}$.

5.3 Ionizing Photon Production Efficiency

We are able to derive the ionizing photon production efficiency from SED fitting by computing for each template used in this work the Lyman continuum photon production rate $N_{\mathrm{LyC}}$ and the observed monochromatic UV luminosity $L_{\nu}$. The ionizing photon production is then $\xi_{\mathrm{ion}}=N_{\mathrm{LyC}}/L_{\nu}$. For our final sample, we find $\log(\xi_{\mathrm{ion}}/\mathrm{erg}^{-1}\hskip 2.84526pt\mathrm{Hz})=26.29^{+0.40}_{-0.38}$ assuming a Calzetti dust attenuation curve and $\log(\xi_{\mathrm{ion}}/\mathrm{erg}^{-1}\hskip 2.84526pt\mathrm{Hz})=26.07^{+0.27}_{-0.30}$ assuming an SMC dust attenuation curve. In most lower redshift studies (e.g., Shivaei et al., 2018), $\xi_{\mathrm{ion}}$ is inferred from a dust corrected Hydrogen line (e.g., H$\alpha$) for which the flux depends mostly on the Lyman continuum photon production rate (e.g., Storey & Hummer, 1995). Given the large number of unconstrained parameters going into our analysis (e.g., attenuation curve), we consider that our result set a lower limit to the average $z\sim 8$ ionizing photon production efficiency with $\log(\xi_{\mathrm{ion}}/\mathrm{erg}^{-1}\hskip 2.84526pt\mathrm{Hz})\geq 25.77$.

We note that in some studies (e.g., Bouwens et al., 2016), $\xi_{\mathrm{ion}}$ is an intrinsic quantity since it is derived by using a dust-corrected UV luminosity, while in our work we derive an observed $\xi_{\mathrm{ion}}$ value since we do not correct the observed UV luminosity for dust. Due to this dust correction, the intrinsic $\xi_{\mathrm{ion}}$ sets a lower limit for the observed $\xi_{\mathrm{ion}}$. However, thanks to the small dust attenuation that we find for our sample, the difference between intrinsic and observed $\xi_{\mathrm{ion}}$ should be small.

The constraints that we obtain on the ionizing photon production efficiency at $z\sim 8$ are consistent with results obtained for the bluest (i.e., least dust attenuated) LBGs and Lyman-$\alpha$ Emitters at $z\sim 2$ (Shivaei et al., 2018; Sobral et al., 2018; Tang et al., 2018) as well as results obtained for low-redshift compact star-forming galaxies (Izotov et al., 2017), and for LBGs at $z\sim 4-5$ (Bouwens et al., 2016; Lam et al., 2019; Ceverino et al., 2019). The observed ionizing photon production efficiency that we find is also consistent with the observed value found for $z\sim 0.3$ Lyman continuum emitters ($\xi_{\mathrm{ion}}=25.6-26$, Schaerer et al., 2016).

This is the first time that the ionizing photon production efficiency is estimated for a significant sample of galaxies in the reionization era and this value is higher than the canonical value $\log(\xi_{\mathrm{ion}}/\mathrm{erg}^{-1}\hskip 2.84526pt\mathrm{Hz})=25.2-25.3$ by a factor $\geq 3$. This higher value of $\xi_{\mathrm{ion}}$ translates into a lower value of the Lyman continuum escape fraction required in a scenario where star-forming galaxies are driving cosmic reionization (Bouwens et al., 2016; Shivaei et al., 2018; Chevallard et al., 2018b; Matthee et al., 2017b, a; Lam et al., 2019).

5.4 Evolution of the EW(H$\beta$+[O III]$\lambda\lambda 4959,5007$) with Redshift

We show the EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$) distribution for our $z\sim 8$ sample in Fig. 4 along with the EW distribution at $z=3.24$ (Khostovan et al., 2016) and the one at $z\sim 6.8$ for the extreme emitter sample of Smit et al. (2015). Our distribution is consistent with the latter, given that the sample of Smit et al. (2015) only included sources with the largest EW. We find a median $\mathrm{EW}([\textrm{O}~{}\textsc{iii}]\lambda\lambda 4959,5007+{\rm H}\beta)=649^{+92}_{-49}$Å, consistent with the value of 670${}^{+260}_{-170}$Å from Labbé et al. (2013).

Comparing our EW distribution with the Sloan Digital Sky Survey (SDSS, Abolfathi et al., 2018), we find that only $0.23\pm 0.01\%$ galaxies exhibit such strong emission lines ($\mathrm{EW}([\textrm{O}~{}\textsc{iii}]\lambda\lambda 4959,5007+{\rm H}\beta)\geq 300$Å) in the entire SDSS sample. We also compare the $z\sim 8$ distribution with two distributions drawn from SDSS. The SDSS samples were mass- or SFR-matched by randomly picking 50 SDSS galaxies within 0.05 dex and 0.1 dex in terms of stellar mass and SFR, respectively, for each galaxy in our sample. The SDSS sample selected through SFR-matching exhibits only a small overlap with the EWs derived in our work ($1.6\pm 0.2\%$), while the one matched by stellar mass leads to a non-negligible fraction of galaxies with EW as high as found in our $z\sim 8$ sample ($9.8\pm 0.4\%$). Nevertheless, it is clear that the emission lines of $z\sim 8$ galaxies are much more extreme than local galaxies of similar mass or SFR, although specific local and low-z population such as Blue Compact Dwarf galaxies exhibit similar properties (e.g., Izotov et al., 2011; Cardamone et al., 2009; Yang et al., 2017; Rigby et al., 2015; Senchyna et al., 2017).

We have also attempted to match our sample in terms of sSFR but the number of SDSS galaxies to exhibit similarly high sSFR as in our sample (median $\mathrm{sSFR}=63^{+188}_{-55}$Gyr*-1*) is extremely small ($<0.35\%$), such that no representative sSFR-matched sample could be constructed. However, the small sample of galaxies with such high sSFR does indeed exhibit EWs as large as derived in our work. This illustrates that finding a significant sample of local galaxies with properties (stellar mass, SFR, sSFR, and EW([O iii]$\lambda\lambda 4959,5007$+H$\beta$) similar to $z\sim 8$ galaxy properties is a difficult task.

From $z=3.24$ to $z\sim 8$, there is a clear evolution of the median of the EW distribution. However, the EWs have to be compared for a given stellar mass range (e.g., $9.5<\log(\mathrm{M}_{\star}/\mathrm{M}_{\odot})<10.0$, Khostovan et al., 2016). The derived stellar mass for our sample is significantly lower than the $z\sim 3.24$ sample, with only four galaxies (5%) with $\log(\mathrm{M}_{\star}/\mathrm{M}_{\odot})\geq 9.5$. A comparison of the EW properties as a function of stellar mass is shown in Fig. 5. While individual error bars are relatively large ($\sim\pm 0.4$ dex for both parameters), our sample is consistent with the $z=3.24$ and $2.23$ EW–M*⋆* relations. Furthermore, using the median stellar mass of our sample with the $z=2.23$ relation from Khostovan et al. (2016), we predict a median equivalent width $\mathrm{EW}=712^{+70}_{-62}$Å for our sample, a value consistent with our derivation based on SED fitting.

6 The H$\beta$+[O iii] Luminosity Function at $z\sim 8$

Based on the photometric estimates of the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ emission line strengths of all the $z\sim 8$ galaxies in the GREATS sample, we have the opportunity for a first derivation of the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ luminosity function at these redshifts, which we describe in the next section.

6.1 Derivation of the Emission Line Luminosity Function

Our approach is based on converting the UV LF to an emission line LF using the relation between the observed UV luminosity, $L_{UV}$, and the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ line luminosity, $L_{OIII+H\beta}$, at $z\sim 8$. This approach is analogous to the one used in the derivation of stellar mass functions at very high redshift (e.g., in González et al., 2010; Song et al., 2016) or the star-formation rate function (Smit et al., 2012, 2016; Mashian et al., 2015).

The relation between the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ luminosity and the observed UV luminosity in our sample is calibrated in Fig. 6. To increase the range of UV luminosities probed in our work, we add to our sample with $S/N(3.6\mu\mathrm{m}\lor 4.5\mu\mathrm{m})\geq 3$ galaxies with lower $S/N$. We apply to these galaxies the same procedure as the rest of the sample and so we obtain the complete probability distribution function for all the parameters, including the [O iii]+H$\beta$ luminosity. While uncertainties remain relatively large for individual galaxies ($\sim 0.18$ dex and $0.19-0.24$ dex for the UV and [O iii]$\lambda\lambda 4959,5007$+H$\beta$ luminosities, respectively), we find that the UV luminosity and the [O iii]+H$\beta$ luminosity are well correlated (Spearman rank correlation coefficient $\rho=0.56$, standard deviation from null hypothesis $\sigma=7.1$).

We use a MCMC method to fit the relation with three parameters, a slope and an intercept, plus an intrinsic (Gaussian) dispersion around the relation, $\sigma_{int}$. This results in:

[TABLE]

Together, with an intrinsic dispersion of $\sigma_{int}=0.35$ dex around the median fit. The corresponding 68-percentile contours are also shown in Fig. 6.

As a comparison, we also use the publicly available mock catalog JAdes extraGalactic Ultradeep Artificial Realizations (JAGUAR, Williams et al., 2018) to derive the relation between UV and [O iii]+H$\beta$ luminosity of simulated $z\sim 8$ galaxies. The JAGUAR mock catalog has been produced by matching luminosity and stellar mass functions as well as the relation between the stellar mass and UV luminosity, mostly at $z\leq 4$. The galaxy properties are then extrapolated up to $z\sim 15$. The JAGUAR catalog provides emission line fluxes and EWs for the main lines based on modeling with the BEAGLE code (Chevallard & Charlot, 2016; Chevallard et al., 2018a). We identify all galaxies from the fiducial JAGUAR mock in the redshift range $7.11<z<9.05$ and we randomly select 1000 of them to match the absolute UV magnitude distribution of our sample, and then fit the UV-[O iii]+H$\beta$ luminosity data. The result is shown in red in Fig. 6. Similarly to our sample, the $z\sim 8$ galaxies from the JAGUAR catalog exhibit a tight relation between UV and [O iii]+H$\beta$ luminosity (Spearman rank correlation coefficient $\rho=0.73$, standard deviation from null hypothesis $\sigma>40$). However, the mock galaxies exhibit a significantly lower [O iii]+H$\beta$ luminosity ($\sim 0.5$ dex) at a given $L_{\mathrm{UV}}$ compared to the relation of our galaxies. The detailed reason for this discrepancy relative to the JAGUAR mock is unclear, but one possible reason is differences in the median physical properties. For instance, while the mock galaxies exhibit (3.6-4.5)µm color similar to the ones from our sample at a given UV luminosity, the average F125W-3.6µm color in JAGUAR is smaller by $\sim 0.3$ magnitude compared to the observed F125W-3.6µm color in our sample. This means that while (3.6-4.5)µm color and EW([O iii]+H$\beta$) are on average similar between JAGUAR and our sample, the absolute [O iii]+H$\beta$ line luminosity scales with the 3.6µm flux which is larger in our sample compared to the JAGUAR mock catalog. Furthermore, JAGUAR models a small field comparatively to our data, therefore the overlap in UV luminosity is small.

Using the relation between the UV and [O iii]+H$\beta$ luminosity of Eq. 6.1, we can now derive the [O iii]+H$\beta$ LF based on the known $z\sim 8$ UV LF (Bouwens et al., 2015). In order to properly compute errorbars, we adopt a Markov Chain Monte Carlo approach (see Sharma, 2017, for a review). In particular, we sample 100’000 points from the UV LF and convert their corresponding $L_{UV}$ values to $L_{OIII+H\beta}$, based on the relation derived above including the appropriate dispersion $\sigma_{int}$. Finally, we fit a Schechter function to the resulting $L_{OIII+H\beta}$ values, keeping the three quantities $\Phi^{\star}$, $\mathrm{L}^{\star}$, and $\alpha$ as free parameters. We repeat this procedure 10’000 times, and vary the input Schechter function parameters of the UV LF according to their appropriate covariance matrix, which was derived from the contour plots shown in Bouwens et al. (2015). This procedure results in 10’000 line LFs, from which we compute the mean and standard deviation for all three Schechter function parameters of the line LF.

The final result is shown in Fig. 7, with Schechter function parameters of the line LF of $\log(\mathrm{L}^{\star}/\mathrm{erg}\hskip 2.84526pt\mathrm{s}^{-1})=43.45^{+0.21}_{-0.19}$, $\log(\Phi^{\star}/\mathrm{Mpc}^{-3})=-4.05^{+0.40}_{-0.44}$, and $\alpha=-2.22^{+0.28}_{-0.32}$. The points with errorbars in Fig. 7 correspond to the observed UV LF measurements from Bouwens et al. (2015), which were converted to the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ LF using the same approach as described above. They clearly agree very well with the mean Schechter function derivation.

6.2 The Evolution of the $[OIII]+H\beta$ Line Luminosity Function to $z\sim 8$

The [O iii]$\lambda\lambda 4959,5007$+H$\beta$ LF has previously been derived up to $z\sim 3$ by several authors (Hippelein et al., 2003; Ly et al., 2007; Colbert et al., 2013; Pirzkal et al., 2013; Khostovan et al., 2015). By adding our estimate at $z\sim 8$, we can thus study its evolution across more than 13 Gyr. Fig. 7 shows such a comparison of the $L_{OIII+H\beta}$ LFs derived at different redshift, from $z\sim 0.5$ to $z\sim 8$. Interestingly, the line LF is found to be higher at all luminosities at $z\sim 8$ compared to $z\sim 3$. This is in stark contrast to the evolution of the UV LF, which peaks at $z\sim 2-3$, but then steadily declines to higher (or lower) redshift. We show in Fig. 8 the evolution with redshift of L⋆ and $\Phi^{\star}$ and by extrapolating the observed trends at lower redshift to $z\sim 8$, especially the Khostovan et al. (2015) results, the $z\sim 8$ values are in remarkable agreement with expectations.

The evolution of the [O iii]+H$\beta$ LF from $z\sim 3$ to $z\sim 8$ can be explained by an evolution of the relation between L(UV) and the stellar mass with redshift. It is known that high-redshift galaxies have their stellar mass related to their UV luminosity (Stark et al., 2009; González et al., 2011; Duncan et al., 2014; Grazian et al., 2015; Song et al., 2016; Stefanon et al., 2017) and while at $0<z\leq 4$ the slope of the $M_{\mathrm{UV}}-\mathrm{M}_{\star}$ relation is not evolving, the intercept evolves in such way that at a given stellar mass the UV luminosity decreases with increasing redshift (Williams et al., 2018). The evolution of this relation is more uncertain at $z>4$ because of the uncertainties on the stellar mass estimation due to the nebular emission contamination (e.g., De Barros et al., 2014), but assuming a similar evolution at $z>4$, combining the evolution of the UV luminosity at a given stellar mass and the relation observed between nebular emission line EW and stellar mass (Fig. 5; Fumagalli et al., 2012; Sobral et al., 2014; Khostovan et al., 2016), we expect that with increasing redshift, at a given UV luminosity, the stellar mass decreases and the EW([O iii]+H$\beta$) increases accordingly. This means that at a given UV luminosity the [O iii]+H$\beta$ luminosity is increasing with increasing redshift. This is the expected trend but uncertainties on the stellar mass estimation (Fig. 5) precludes any further quantification of the UV luminosity vs. stellar mass relation evolution from $z\sim 3$ to $z\sim 8$.

To test this explanation, we derive the relation between L([O iii]+H$\beta$) and L(UV) at $z\sim 3$ through abundance matching. In particular, we use the $z\sim 3$ UV LF from Reddy & Steidel (2009) to derive the cumulative number density of galaxies at a given UV luminosity and match it to the corresponding cumulative number density at a given [O iii]+H$\beta$ luminosity based on the $z=3.2$ [O iii]+H$\beta$ LF from Khostovan et al. (2015). We show the resulting relation in Fig. 6 with a blue line. There is a clear evolution of the L(UV) vs. L([O iii]+H$\beta$) relation from $z\sim 3$ to $z\sim 8$, galaxies being brighter in [O iii]+H$\beta$ at any given UV luminosity explored in this work and an increasing difference between the $z\sim 3$ and $z\sim 8$ relation with decreasing UV luminosity. This finding supports our explanation for the [O iii]+H$\beta$ LF evolution.

6.3 Predictions of JWST Number Counts

One of the most awaited capabilities of the upcoming $JWST$ is its unprecedented sensitivity for spectroscopy at $>2\micron$. In particular, $JWST$ will for the first time provide spectroscopic access to the rest-frame optical emission lines of very high redshift galaxies, including the [O iii]$\lambda\lambda 4959,5007$ and H$\beta$ lines at $z\sim 8$ that we constrained through photometry here. In order to decide on the area and depth for the most efficient spectroscopic surveys with $JWST$, an estimate of the [O iii]$\lambda\lambda 4959,5007$+H$\beta$ LF as derived above is of critical importance.

Of particular interest for such $JWST$ predictions is the NIRSpec instrument (Bagnasco et al., 2007), which will be the workhorse NIR spectrograph. With four quadrants, each of which covers 2.3 arcmin2, NIRSpec spans an effective area of $\sim 9.2$ arcmin2. Its sensitivity is exquisite. In only 1hr, NIRSpec will reach $6\sigma$ detections for emission lines at $\sim$4.5$\micron$ and fluxes of $\sim$6$\times 10^{-19}$ erg s*-1cm-2* at $R=1000$, or $\sim$8$\times 10^{-19}$ erg s*-1cm-2* at $R=100$. These numbers were derived with the latest JWST/ETC, when integrating over the full extent of the lines (which were assumed to have an intrinsic width of 150 km s*-1*).

In Fig. 9, we plot the cumulative surface density of galaxies at $z=7.5-8.5$ as a function of emission line luminosities based on the line LF which we derived in the previous section. In particular, we show the total combined [O iii]$\lambda\lambda 4959,5007$+H$\beta$ luminosities (as would be seen, e.g., in $R=100$ low-resolution spectroscopy with NIRSpec), and we also split the luminosities into the three different lines, $H\beta$, [OIII]4959, and [OIII]5007. For the latter step, we employ the median line ratios as found in the SED fitting in Section 4, with H$\beta$/(H$\beta$+[O iii]$\lambda\lambda 4959,5007$)=0.21, [O iii]$\lambda 4959$/(H$\beta$+[O iii]$\lambda\lambda 4959,5007$)=0.20, and [O iii]$\lambda 5007$/(H$\beta$+[O iii]$\lambda\lambda 4959,5007$)=0.59.

As can be seen, at the 1hr sensitivity limits of NIRSpec with $R=100$, we can expect a cumulative surface density of $z=7.5-8.5$ galaxies of 101.2±0.3 arcmin*-2* which have blended [OIII]$+H\beta$ lines that are bright enough to be significantly detected. This means that a single NIRSpec mask (with effective area 9.2 arcmin2) would, in principle, be able to target on average $\sim 150$ galaxies (where the $1\sigma$ uncertainties range from 80 to 300 galaxies). The fixed grid of the NIRSpec slitlet masks will reduce this number somewhat. Unfortunately, however, there is an additional limitation of early spectroscopic surveys. The depth reached in terms of [O iii]+H$\beta$ luminosity in 1hr for NIRSpec corresponds to an observed UV magnitude of $m_{\mathrm{UV}}=29.9$ ($\sim 29-31$ when accounting for the uncertainties in the UV-[O iii]+H$\beta$ luminosity relation, Eq. 6.1). This implies that the average surface density of current $z\sim 8$ galaxy samples in the prime extragalactic legacy fields such as CANDELS is significantly lower than the number above (Bouwens et al., 2015). Therefore, early JWST spectroscopic surveys, which are based on the selection of targets from current $HST$ datasets, cannot be maximally efficient for a targeted $z\sim 8$ galaxy survey. The best strategy will be to perform deep pre-imaging with JWST to identify targets at $z\sim 8$.

Nevertheless, our calculation shows that if significantly deep imaging data are available to select targets from, a single NIRSpec mask with $R=100$ could be filled with just $z\sim 8$ emission line sources for which a 1hr observation can measure a secure redshift. This will result in revolutionary insights of the large scale structure in the heart of the reionization epoch.

Of course, in order to study the physics, more than a simple redshift measurement is required. In particular, the [OIII]$+H\beta$ lines need to be split with observations at $R=1000$ or higher. For such surveys, the corresponding number of expected galaxies in 1hr observations are only 24 galaxies for $H\beta$ lines, and 97 galaxies with 6$\sigma$ [OIII]5007 line detections, per NIRSpec mask.

As a final remark, we compare our predictions with the ones from the JAGUAR mock catalog. Given the lower [O iii]+H$\beta$ luminosities at a given L(UV) compared to our observed sample, it is clear that the mocks will significantly underpredict the number of observed sources at a given line luminosity. When computing the line LF from the JAGUAR catalogs, we find a lower normalization by a factor $\sim 8\times$ compared to our observed LF. Hence, the JAGUAR mock catalogs should underpredict the number of rest-frame optical emission lines that can be detected at $z\sim 8$ with $JWST$/NIRSpec in the future by the same factor.

7 Conclusions

We have presented a detailed analysis of a $z\sim 8$ galaxy sample with some of the deepest available Spitzer observations, from the GREATS survey. The sample has been culled through photometric redshifts to ensure that the selected galaxies are at $z\geq 7.11$, where the IRAC 3.6µm$-$IRAC 4.5µm colors put strong constraints on the [O iii]+H$\beta$ equivalent width. We built and used a photoionization grid with a large parameter space covering a variety of stellar metallicity and ISM conditions, using the BPASS models as stellar emission inputs. We feed the resulting SEDs that include stellar and nebular emission (continuum and lines) to our SED fitting code to derive $z\sim 8$ galaxy properties. Accounting for the photometric and model uncertainties, we have specifically derived the [O iii]+H$\beta$ luminosities, allowing us to derive the corresponding [O iii]+H$\beta$ luminosity function and make predictions for JWST observations.

In summary, we find the following.

Our subsample with $S/N(3.6\mu\mathrm{m}\lor 4.5\mu\mathrm{m})\geq 3$ has the following average properties: $\log(\mathrm{M}_{\star}/\mathrm{M}_{\odot})=8.62^{+0.43}_{-0.39}$, $\log(\mathrm{age}/\mathrm{yr})=7.2^{+0.9}_{-0.6}$, $\log(\mathrm{SFR}/\mathrm{M}_{\odot}{\rm yr}^{-1})=1.26^{+0.42}_{-0.30}$, $A_{V}=0.4\pm 0.2$, and $\mathrm{sSFR}=63^{+188}_{-55}\mathrm{Gyr}^{-1}$. 2. 2.

To reproduce the observed IRAC color of this subsample, which is strongly affected by [O iii]$\lambda\lambda 4959,5007$+H$\beta$ emission, the two main parameters driving the EW are the stellar metallicity and the ionization parameter, and they have the following values $Z_{\star}=0.004_{-0.002}^{+0.004}$ and $\log\mathrm{U}=-3.0\pm 1.0$. 3. 3.

We are able to put constraints on the median ionizing photon production efficiency with $\log(\xi_{\mathrm{ion}}/\mathrm{erg}^{-1}\hskip 2.84526pt\mathrm{Hz})\geq 25.77$. This latter value is $>3$ times higher than the canonical value, implying that these galaxies have a higher ionizing output than typically assumed and can thus more easily reionize the universe. 4. 4.

According to our SED fitting which matches the observed IRAC colors, we find a median rest-frame equivalent width $\mathrm{EW}([\textrm{O}~{}\textsc{iii}]\lambda\lambda 4959,5007+{\rm H}\beta)=649^{+92}_{-49}$Å (Fig. 5). 5. 5.

We find a relatively tight relation between [O iii]+H$\beta$ and UV luminosity (Fig. 6), allowing us to derive for the first time the [O iii]+H$\beta$ LF at $z\sim 8$ based on the $z\sim 8$ UV LF. We find that, in contrast with the evolution of the UV LF from $z\sim 3$ to $z\sim 8$, the $z\sim 8$ [O iii]+H$\beta$ LF is higher at all luminosities than at $z\sim 3$. This is due to the increasing [O iii]+H$\beta$ luminosity at a given UV luminosity with increasing redshift. 6. 6.

Finally, we use the derived [O iii]+H$\beta$ LF to predict JWST number counts. A single NIRSpec pointing would contain $\sim 150$ galaxies at $z=7.5-8.5$, for which the [O iii]+H$\beta$ emission could be detected in only 1hr. However, the current average surface density of $z\sim 8$ galaxies in extragalactic legacy fields is significantly lower that this number. Therefore, to maximize the efficiency of JWST  deep pre-imaging to $m_{\mathrm{UV}}\sim 30$ mag will be required.

In this work, we have used the deepest Spitzer data available on large areas. We have accounted for observational uncertainties on the photometry and we have used a grid of photoionization models with a large parameter space. While we have attempted to minimize the number of assumptions going into our analysis, many modeling uncertainties are still present: the ingredients in the stellar population synthesis models (stellar atmospheres, binaries, rotation), the IMF, the possible presence of multiple stellar populations, the dust attenuation curve, the ratio between nebular and stellar attenuation, interstellar abundances, and depletion factor of metals on to dust grains. Only the unprecedented abilities of JWST will allow to alleviate some of these uncertainties.

Acknowledgements

We thank the anonymous referee who helped improve this manuscript. The work of SDB has been partially supported by a Flexibility Grant from the Swiss National Science Foundation and by a MERAC Funding and Travel Award from the Swiss Society for Astrophysics and Astronomy. V. G. was supported by CONICYT/FONDECYT initiation grant number 11160832.

This work made use of v2.1 of the Binary Population and Spectra Synthesis (BPASS) models as last described in Eldridge et al. (2017). Calculations were performed with version 17.00 of Cloudy  last described by Ferland et al. (2017). This work also made use of Astropy, a community-developed core Python package for Astronomy (The Astropy Collaboration et al., 2018), as well as the pymc3 library (Salvatier et al., 2016).

This paper made use of public catalogs derived from data taken by the Sloan Digital Sky Survey IV. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.