Elevated ionizing photon production efficiency in faint high-equivalent-width Lyman-α emitters

ABSTRACT

While low-luminosity galaxies dominate number counts at all redshifts, their contribution to cosmic reionization is poorly understood due to a lack of knowledge of their physical properties. We isolate a sample of 35 z ≈ 4–5 continuum-faint Lyman-α emitters from deep VLT/MUSE spectroscopy and directly measure their H α emission using stacked Spitzer/IRAC Ch. 1 photometry. Based on Hubble Space Telescope imaging, we determine that the average UV continuum magnitude is fainter than −16 (≈ 0.01 L^⋆), implying a median Lyman-α equivalent width of 259 Å. By combining the H α measurement with the UV magnitude, we determine the ionizing photon production efficiency, ξ_ion, a first for such faint galaxies. The measurement of log₁₀ (ξ_ion [Hz erg⁻¹]) = 26.28 (⁠|$^{+0.28}_{-0.40}$|⁠) is in excess of literature measurements of both continuum- and emission line-selected samples, implying a more efficient production of ionizing photons in these lower luminosity, Lyman-α-selected systems. We conclude that this elevated efficiency can be explained by stellar populations with metallicities between 4 × 10⁻⁴ and 0.008, with light-weighted ages less than 3 Myr.

Galaxies: evolution, Galaxies: high-redshift, Galaxies: ISM

1 INTRODUCTION

Although recent observations have provided ever more certainty about the timing and duration of the last significant phase transition in the Universe, cosmic reionization (e.g Bañados et al. 2018; Planck Collaboration 2018), much remains to be understood about the source(s) of the photons that caused the change. Evidence is mounting that star-forming galaxies could have produced enough ionizing photons to drive reionization at z > 6 but typically only under the assumption that galaxies with UV magnitudes much fainter than the characteristic luminosity (L^⋆) dominate the total number counts and that these galaxies have similar physical properties to brighter systems (e.g. Yan & Windhorst 2004; Bouwens et al. 2012; Finkelstein et al. 2012; Robertson et al. 2013).

While observed UV luminosity functions indicate that faint galaxies are numerous at all redshifts (e.g. Atek et al. 2018), few direct observational (spectroscopic) constraints on the efficiency of their ionizing photon production exist. This difficulty is primarily due to the observability of spectral features, particularly those in the rest frame optical: ground-based spectroscopy can only cover features such as H α in the near-infrared (IR) until z ≈ 2.8. H α, in particular, is crucial to understand the contribution of galaxies to reionization as it is directly related to the intrinsic production rate of ionizing photons (compared to the resonantly scattered Lyman-α). The ratio of the H α flux to the flux of non-ionizing (UV) photons, when combined with the escape fraction of ionizing photons and the number density of galaxies, can be used to determine the total production rate of ionizing photons in the early Universe.

Although detecting H α is currently not possible with traditional spectroscopy for most systems at z > 2.8 (and detections of H β are often infeasible due to its faintness in low-luminosity galaxies), photometric techniques have been developed to measure H α and other rest frame-optical emission lines in longer wavelength imaging data. For example, the existence of bright, high-equivalent-width (EW) optical emission lines in z ≳ 4 galaxies is inferred by measuring strong excesses in broad-band Spitzer/IRAC Ch.1 (3.6 μm) and/or Ch. 2 (4.5 μm) photometry (e.g. Shim et al. 2011; González et al. 2012; Labbé et al. 2013; Smit et al. 2014; Rasappu et al. 2016; Roberts-Borsani et al. 2016). These studies have demonstrated that the typical rest frame EWs of [O iii] and H α in ∼L^⋆ galaxies at high z often exceed 300 Å. Moreover, the highest EW sources at z > 7 have ionizing photon production efficiencies that are a factor of 2.5 higher than the ‘canonical’ value, implying a diversity in the ionization properties of the full galaxy population and that strong line emitters could be important contributors to cosmic reionization (Matthee et al. 2017b; Stark et al. 2017).

The vast majority of the total galaxy population, namely those at sub-L^⋆ UV luminosities, are much more difficult to understand at these redshifts. The low spatial resolution (FWHM > 1.5 arcsec) and broad filter widths (>6000 Å) make shallow IRAC spectro-photometry for individual objects useful only when they are isolated and have relatively bright emission lines. Photometric stacking, though, can be used in order to detect fainter emission lines so long as there are no nearby contaminating sources, or those sources can be accurately modelled. As shown in Lam et al. (2019), this work can be extended to faint (M_UV > −18; 0.05 L^⋆ at z = 4) galaxies when using MUSE spectroscopy and ultra-deep Spitzer/IRAC imaging with new techniques for modelling contamination from nearby sources (Labbé et al. 2015).

Moving to even fainter sources is necessary to fully understand the ‘budget’ of ionizing photons in the early Universe (Finkelstein et al. 2019). Sources fainter than M_UV ≈ −17 do not appear in the deepest broad-band imaging taken with Hubble unless they have been gravitationally lensed. Yet, an abundant population of sources with M_UV ≈ −15 (0.01 L^⋆) at z ≈ 2.9−6.7 has recently been discovered spectroscopically via Lyman-α emission in un-targeted surveys with the MUSE spectrograph (Bacon et al. 2017; Maseda et al. 2018). Their bright Lyman-α emission and faint UV continuum (detected only in stacks) implies that they have extreme Lyman-α EWs, in excess of 150 Å. Their number density is consistent with a simple extrapolation from the general population of star-forming Lyman-α-emitters at these redshifts. Their high EWs can be produced only in stellar populations with ages less than 10 Myr and metallicities less than a few per cent Z_⊙ (e.g. Malhotra & Rhoads 2002,Raiter, Schaerer & Fosbury 2010; Hashimoto et al. 2017a, cf. Neufeld 1991), physical conditions which are also conducive to efficiently producing ionizing photons.

Here we aim to study the ionizing photon production efficiency in these high-EW Lyman-α emitters (LAEs) by using stacked Spitzer/IRAC photometry to measure the H α emission. We can put these results into context by comparing to the numerous studies at similar redshifts that have probed more luminous galaxies: LAEs (e.g. Harikane et al. 2018; Lam et al. 2019), H α emitters (HAEs; Matthee et al. 2017a), and continuum dropouts (e.g. Bouwens et al. 2016). Our values will also be compared to larger samples of ‘typical’ ≈ 0.3−3 L^⋆ star-forming galaxies at z ≈ 2 (e.g. Shivaei et al. 2018) and theoretical stellar population models for the evolution of ξ_ion with physical properties.

This article is organized as follows. In Section 2, we describe the MUSE and IRAC data sets used for this study. In Section 3, we discuss how the data are used to determine the H α luminosity and hence ξ_ion, with a discussion of these results given in Section 4. Finally, in Section 5, we summarize the primary results and give an outlook for future studies. We adopt a flat ΛCDM cosmology (Ω_m = 0.3, |$\Omega _\Lambda =0.7$|⁠, and H₀ = 70 km s⁻¹ Mpc⁻¹) and AB magnitudes (Oke 1974) throughout.

2 DATA

Using data from the MUSE spectrograph (Bacon et al. 2010) on the Very Large Telescope, we select high-EW LAEs from the MUSE UDF survey (Bacon et al. 2017; Inami et al. 2017). This survey covers approximately 9 arcmin² to a depth of 10–30 h at optical wavelengths (4750–9300 Å). The un-targeted nature of MUSE spectroscopy means that objects do not need to be pre-selected based on their continuum properties, and we can search the data cubes for spectral features such as emission lines without requiring the detection of a continuum counterpart. From these catalogues, we consider only LAEs in the range 3.829 < z < 4.955 with confident redshifts (i.e. CONFID of 2 or greater), where H α emission lies within the IRAC Ch. 1 bandpass.

In order to select a sample of high-EW LAEs, we must select not only on the Lyman-α flux but also on the UV continuum level. Using the aperture flux measurements from Maseda et al. (2018), based on the combined 11-band Illingworth et al. (2013) and Oesch et al. (2018) reductions, we select only MUSE LAEs from the Inami et al. (2017) catalogue that have no >3-σ HST detections in bands redwards of Lyman-α. To wit, objects detected only in the HST band that contains Lyman-α are still included in the sample. This results in a sample of 41 objects. This selection differs from that presented in Lam et al. (2019), which is also based on Lyman-α-emitters from the MUSE UDF survey, as they require continuum detections for individual objects. The threshold of 3-σ compared with 5-σ as in Maseda et al. (2018) is to restrict ourselves to only the highest-EW sources by limiting the acceptable detection limit in the UV continuum (the average F850LP 3-σ magnitude limit in the field is approximately 30.0; Illingworth et al. 2013). We note that the qualitative conclusions of this work are independent of this choice.

A key addition to the HST photometry is near-IR observations with Spitzer/IRAC. Here we utilize the combination of data from the ‘GOODS Re-ionization Era wide-Area Treasury from Spitzer’ (GREATS, PI: I. Labbé; M. Stefanon, in preparation) as well as all previous IRAC observations taken over the larger area around the UDF (Labbé et al. 2015; Lam et al. 2019). The average exposure time of the IRAC Ch. 1 imaging is 167 h (26.6 magnitude using the photometric procedure of Labbé et al. 2015, 5-σ), with some regions in the UDF receiving 278 h of coverage. In IRAC Ch. 2, the exposure time is similar (average 139 h, deepest 264 h, 26.5 magnitude).

From our candidate sample of 41, we manually remove six objects with severely contaminated IRAC photometry and those sources with significant residual flux in the cleaned images (see Section 2.1). This leaves a final sample of 35 high-EW LAEs with a median redshift of z = 4.52, for which HST/ACS F775W cutouts are shown in Fig. 1. This represents 19 per cent of all LAEs in this redshift range from the MUSE UDF catalogues (Inami et al. 2017, and Bacon et al. in preparation). Based on the UV continuum slope measured for these galaxies described in Section 3.4 and the median Lyman-α flux from the point spread function (PSF)-weighted MUSE spectra (which does not include extended emission; Inami et al. 2017), the average rest frame Lyman-α EW of the sample is 259 Å.¹

Figure 1.

HST/ACS F775W cutouts showing the z ≈ 4–5 LAEs used in this study. Purple cross hairs denote objects that are not detected in any HST photometric band at a significance greater than 3-σ (cf. Maseda et al. 2018), while the blue cross hairs denote an object that is detected only in the HST photometric band that contains Lyman-α emission (i.e. F775W and/or F606W). Cutouts are 5 arcsec on a side, and the length of the arms of the cross hairs is 1 arcsec. The dashed black circle in the upper-left panel shows the FWHM of the IRAC Ch. 1 PSF (1.6 arcsec diameter). A colour version is available online.

2.1 De-blended IRAC photometry

The large spatial PSF of IRAC (FWHM ≈ 2 arcsec) presents a challenge when dealing with compact sources in crowded fields. This is especially true in deep imaging data, where the wings of the PSF from neighbouring sources often overlap. To correct for this effect, a higher spatial resolution image can be combined with a model PSF in order to de-blend the photometry in a crowded field. In this case, we use HST F850LP imaging as the high-resolution prior image and follow the procedure of Labbé et al. (2015) and Lam et al. (2019). We create a model IRAC flux distribution for all continuum-detected sources in the field, which we then subtract from the IRAC data. We are therefore left with a residual image that should contain only flux from the primary (continuum-undetected) source

All 12 × 12 arcsec cutouts are visually inspected for defects, which can be caused by, e.g. poorly modelled bright sources or can occur in exceptionally crowded fields. Any object that is deemed to have contaminated photometry due to model residuals is removed from the sample. We note that these contamination effects are independent of the intrinsic properties of the primary source. The resulting ‘clean’ sample has residuals in the range of 0.01 to 0.3 nJy (cf. the mean from the ‘contaminated’ sources of 0.6 nJy), determined via the standard deviation of pixel values within an annular region centred on the source and extending from 1.5 to 3 arcsec.

3 DETERMINATION OF ξ_ION

The three ingredients required to measure ξ_ion are the flux of H α (a proxy for the intrinsic ionizing photon production rate), the flux of the UV, and the escape fraction of ionizing photons:

$$\begin{eqnarray*} \xi _{\textrm{ion}} (\textrm{Hz erg}^{-1})= Q(H^0) / L_{\textrm{UV,int}} , \end{eqnarray*}$$

(1)

where Q(H⁰) is the intrinsic rate of ionizing photons with units of s⁻¹ and the intrinsic UV luminosity, L_{UV, int}, has units of erg s⁻¹ Hz⁻¹. Below we outline the determinations of each of these quantities from our data.

3.1 H α from IRAC stacking

We create a three-dimensional array of the cleaned IRAC cutouts for the LAEs, each of which is centred on the peak of the Lyman-α flux. We calculate the mean of the array at each spatial pixel, incorporating a sigma-clipping procedure with a threshold of 2-σ in order to ensure that individual bright pixels do not dominate the mean. This ‘stacking’ procedure results in a two-dimensional image, representing an average of the input sources. As the IRAC PSF varies across the field due to the combined nature of the GREATS data set (see section 3.2 of Labbé et al. 2015), we also create a three-dimensional array of the local PSFs at the position of each LAE. This stack of PSFs is combined in the same way as the stack of the science data arrays, using the same mask derived from sigma clipping. We fit this PSF to the stacked image to determine the total flux, under the assumption that our sources are unresolved compared to the 2 arcsec IRAC PSF (cf. the HST images). The fit includes uncertainties in the flux determined from stacking the IRAC noise images in the same way as the data. This stacking procedure is used to measure the photometry and create the images shown in Fig. 2.

Figure 2.

Top: Mean-stacked HST, K_s, and IRAC photometry for the 35 high-EW LAEs in the sample. The HST/K_s cutouts are 5 arcsec on a side, while the IRAC cutouts are 12 arcsec. The apertures shown are used for the photometric measurements of M_UV from HST (0.4 arcsec radius). Bottom: Total observed rest frame photometry and sample properties of the high-EW LAEs used in this study, as described in the text. Upper limits (3-σ) are denoted with downward-facing triangles. For illustrative purposes, we also plot an arbitrarily scaled theoretical model for a single stellar population 2 Myr after a burst of star formation, with a metallicity of 0.001 (see Section 4 for details). Lines of hydrogen and helium are included in the model spectrum, as is the average intergalactic medium transmission at this redshift from Inoue et al. (2014).

For the remainder of this work, we deal with the distribution of the output IRAC Ch. 1 and Ch. 2 fluxes using random subsets of the 35 individual LAEs. Namely, we create a series of 10 000 data arrays made up of a random set of 35 LAEs, where we sample the full set with replacement. Each of these arrays is analysed as above to determine flux in Ch. 1 and Ch. 2. This bootstrap procedure allows us to measure the distribution of H α fluxes within the sample while taking into account photometric uncertainties.

We determine the H α fluxes for each of the 10 000 bootstrap iterations by subtracting the Ch. 2 flux from the Ch. 1 flux, assuming that Ch. 2 probes the underlying stellar continuum (see Section 3.1.1). We attribute all of the excess flux in Ch.1 to Hα emission: in the models of Gutkin, Charlot & Bruzual (2016), for metallicities ≲10 per cent Z_⊙ the ratio of [N ii] ([S ii]) to H α never exceeds 0.01 (0.1) over a wide range in physical conditions in the interstellar medium. While a detailed determination of the metallicities of these systems is beyond the scope of this paper and will require additional spectroscopic data, we stress that only models of normal stellar populations with extremely young ages (<3 Myr) and/or low metallicities (Z < 0.008; 0.4 Z_⊙) can reproduce such large Lyman-α EWs (Hashimoto et al. 2017a, b, and Section 4). Similar conclusions are drawn by Trainor et al. (2016) for the faintest z ≈ 2.5 LAEs, where they have spectroscopic access to strong optical emission line ratios.

3.1.1 The continuum level around H α

In order to properly measure the emission line flux from a photometric broad-band measurement, we need to establish the local spectral continuum level. For the aforementioned studies that perform measurements of H α EW at high-z, the continuum level at the position of H α is typically established with photometric detections in the bands bluewards/redwards of H α, such as the K_s-band or IRAC Ch. 2. In our case, the ground-based K_s-band data (from zFOURGE; Straatman et al. 2016) are not deep enough for reliable detections of the continuum flux bluewards of H α. As the IRAC Ch. 2 stacks are deeper and are subject to similar systematics as the Ch. 1 data, we use them to constrain the continuum redwards of H α directly.

We assume a flat continuum slope in f_ν (van der Wel et al. 2011). Any residual contamination, if present in both Ch. 1 and Ch. 2, would not bias our determination of the line flux (and hence ξ_ion) since this is determined via the difference between both measurements. Furthermore, residual contamination would lower the measured H α EW for a spectrum that is flat in f_ν. We also assume no contribution to the Ch. 2 flux from emission lines such as [S iii], which is present for the z < 4.51 subset of LAEs. As shown in Lam et al. (2019), [S iii] λ9069 can be strong in LAEs, with EW values in excess of 100 Å. Constraints on [S iii] do not currently exist for galaxies as faint as the LAEs probed here, but at metallicities below 0.3 Z_⊙ the ratio of [S iii] to H α is predicted to be less than 0.05 (Charlot & Longhetti 2001). Furthermore, any contribution of [S iii] to the IRAC Ch. 2 photometry would mean that we are overestimating the continuum level and hence underestimating the strength of H α (and ξ_ion).

We determine a mean H α EW for the sample of 632 Å. The 68 per cent confidence interval, weighted by the signal to noise of each bootstrap measurement of the H α EW, is between 210 and 1600 Å. The large range in EW measurements is primarily driven by the lower signal to noise in the IRAC Ch. 2 stacks, which is typically a factor of 1.9 lower than that derived in IRAC Ch. 1. In the case of measurements at low-EW (i.e. <210 Å), the typical signal to noise in the EW determination is 2.1.

The mean measured star formation rate (SFR) from the H α luminosity is 1.2 M_⊙ yr⁻¹ (0.4–2 M_⊙ yr⁻¹; 68 per cent) when using the Murphy et al. (2011) conversion. This is larger than the implied UV-based SFR of 0.1 M_⊙ yr⁻¹, likely due to the fact that H α emission probes star formation on shorter time-scales than the UV, namely <10 Myr, which is important given the implied young ages for these systems: see Section 4.

3.2 UV luminosity from HST stacking

Following the same procedure as for the IRAC stacks, for each bootstrap sample we create stacks of the HST images in order to determine M_UV; for more details on this procedure, see Maseda et al. (2018). For the majority of galaxies in our sample, ACS/F775W probes the rest-UV continuum directly. However, for the highest-z objects this band contains flux from Lyman-α: three of the 35 objects have Lyman-α at a position where the throughput of F775W is greater than 33 per cent. The major results of this paper are consistent within the errors (<1-σ) when these higher-z sources are included. Throughout, we assume that F775W probes the UV continuum alone for our sample. Any potential contamination to this flux from Lyman-α emission would mean that we are underestimating ξ_ion. Additionally, an extrapolation assuming a UV continuum slope of −2.5 from the measured F850LP magnitude is consistent with the measured F775W flux (see Section 3.4). Such a continuum slope is expected for faint LAEs (Hashimoto et al. 2017a), and therefore this implies thatF775W is a good tracer of the rest frame UV continuum. We convert the measured M_UV into a luminosity assuming the median redshift for the LAEs in the bootstrap iteration.

3.3 Ionizing photon escape fraction

The escape fraction of ionizing photons relates the observed ionizing photon flux with the intrinsically produced ionizing photon flux:

$$\begin{eqnarray*} f_{\mathrm{esc}}^{\mathrm{ion}} = \frac{Q_{\mathrm{obs}}(H^0)}{Q(H^0)}. \end{eqnarray*}$$

(2)

For brevity, we will refer to |$f_{\mathrm{esc}}^{\mathrm{ion}}$| as f_esc throughout the remainder of the text. This parameter is difficult to measure directly as Q(H⁰) is not a readily observable quantity. Studies such as Ono et al. (2010) suggest using spectral energy distribution (SED) fitting to the broad-band photometry to estimate f_esc indirectly. At z = 6−7, they find f_esc to be consistent with zero and constrain it to be <0.6. Using a similar method, Harikane et al. (2018) fit a relation to f_esc versus EW_Lyα, finding a value consistent with zero for rest frame EWs in excess of 100 Å. We therefore assume f_esc = 0 throughout but note that a higher value of f_esc will result in a larger ξ_ion.

3.4 Dust attenuation

The true amount of dust extinction in these systems is difficult to establish with SED fitting as even in stacked images we observe only the rest-UV portion of the galaxies in addition to the H α emission. Many studies, particularly at these redshifts z ≈ 4–5, instead rely on a measurement of the UV continuum slope β and assume a dust law, as in Meurer, Heckman & Calzetti (1999). We measure a median β = −2.43 from power law fits to the rest-UV photometry from each of the bootstrap iterations, with a 68 per cent confidence interval from −3.00 to −2.04, from a power law fit to the stacked HST photometry for the full stacked sample. We defer a detailed treatment of the β slopes of the full sample of high-EW LAEs to a forthcoming paper, but we would like to highlight that the intrinsic stellar UV continuum slope for systems with the highest Lyman-α EWs must be even bluer than the observed value due to the contribution of nebular continuum emission (Raiter et al. 2010).

For β < −2.23, Meurer et al. (1999) estimate zero dust correction. As our best estimate of β (and 68 per cent of all bootstrap iterations) is below this, we determine that we do not need to include an additional term for dust attenuation when measuring line fluxes. Although our measurement is uncertain due to the intrinsic faintness of the sample, we a priori expect such a correction to be negligible based on results from (UV-brighter) LAEs at similar or higher redshifts (e.g. Stark et al. 2015; Hashimoto et al. 2017a; Harikane et al. 2018) and the result from Trainor et al. (2016), which shows an anticorrelation between nebular reddening and EW_Lyα. Similarly, Tang et al. (2019) measure the dust attenuation via the Balmer decrement for high-EW [O iii] and H α emitters, finding negligible extinction for the highest EW objects. Finally, such a result is also expected when considering the implied low stellar masses and gas-phase metallicities required to power the Lyman-α emission (Garn & Best 2010).

3.5 The sample distribution of ξ_ion

For each bootstrap iteration, we convert the H α flux into a luminosity assuming the median redshift for the LAEs that contributed to the stack. We then convert the luminosity (without a dust correction, as explained above) into the ionizing photon production rate Q(H⁰) assuming Case B recombination at a temperature of 10⁴ K according to:

$$\begin{eqnarray*} \left[\frac{Q(H^0)}{\textrm{s}^{-1}}\right] \times (1-f_{\textrm{esc}}) = \left(\frac{L_{H \alpha }}{\textrm{erg s}^{-1}}\right) \times 7.37 \times 10^{11}, \end{eqnarray*}$$

(3)

as in, e.g. Murphy et al. (2011).

We can therefore calculate ξ_ion according to the usual formula given in equation (1). The weighted mean of the distribution from the bootstrap iterations assuming zero f_esc is log₁₀ (ξ_ion / Hz erg⁻¹) = 26.28, with a 68 per cent confidence interval from 25.89 to 26.56, weighted by the individual uncertainties on each measurement.

In Fig. 3, we show our estimate of ξ_ion versus M_UV, compared to literature samples of line-selected samples (HAEs and LAEs Matthee et al. 2017a; Harikane et al. 2018; Lam et al. 2019) and continuum-selected galaxies (Lyman-break galaxies or LBGs; Bouwens et al. 2016) at similar redshifts, as well as ‘normal’ M_⋆ > 10⁹ M_⊙, L ≈ 0.3–3 L^⋆ star-forming galaxies at z ≈ 2 from Shivaei et al. (2018). We use their value derived from an SMC (Small Magellanic Cloud) extinction curve for consistency with the other studies presented here. The bottom panel shows the relationship between the measured luminosities of H α and the rest frame-UV, which are the two components that go into ξ_ion; any trend in ξ_ion versus M_UV could simply be because the UV luminosity goes into both the ordinate and the abcissa. While the information content in both panels is the same, the axes in the top panel are necessarily correlated.

Figure 3.

ξ_ion versus M_UV for this sample, emission line-selected literature samples of Matthee et al. (2017a), Harikane et al. (2018), and Lam et al. (2019), and the continuum-selected sample of Bouwens et al. (2016). In the case of Shivaei et al. (2018), the dashed line denotes an extrapolation below their faintest bin. The legend describes which result we use from the literature when the authors present multiple values for ξ_ion depending on, e.g. their assumed dust law. The error bars on our data points show the (marginalized) sample distribution based on bootstrapping (see Section 3.1), whereas the error bars on M_UV from the literature sample reflect the bin width used in each study. The bottom panel isolates the two individual components of ξ_ion, as the axes in the top panel are strongly correlated: the quoted uncertainties in ξ_ion for the literature samples do not typically include the width of the M_UV bins. Our high-EW LAEs have higher ionizing photon production efficiency for their continuum (UV) luminosity than the literature samples.

While all of the LAE- and LBG-derived data points are within 2-σ (assuming that the bin widths correspond to a measurement uncertainty on L_UV) of the Shivaei et al. (2018) relation or its extrapolation in either L_Hα or L_UV, the point derived here is more significantly above the relation for L_Hα and L_UV. Therefore, relative to their UV luminosity, the MUSE high-EW LAEs are more efficient at producing ionizing photons than the general population of more luminous LBGs and LAEs at these redshifts.

3.6 The Lyman-α escape fraction

As discussed in Section 3.3, it is difficult to directly measure f_esc in these galaxies, although based on other samples we expect it to be low. We can, however, measure the escape fraction of Lyman-α photons:

$$\begin{eqnarray*} f_{\mathrm{esc}}^{\mathrm{Ly\alpha }} = \frac{L_{\mathrm{Ly\alpha }}^{\mathrm{obs}}}{L_{\mathrm{Ly\alpha }}^{\mathrm{int}}} = \frac{L_{\mathrm{Ly\alpha }}^{\mathrm{obs}}}{8.7 \times L_{\mathrm{H\alpha }}^{\mathrm{int}}}, \end{eqnarray*}$$

(4)

where the superscripts obs and int refer to the observed and intrinsic luminosities, respectively. This equation is valid for Case B recombination with a temperature of 10⁴ K (see Henry et al. 2015; Trainor et al. 2015). Based on our observations, we derive a mean |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$|value of 0.217, with a 68 per cent confidence interval of 0.067–0.333 based on our bootstrap iterations.

4 DISCUSSION

Qualitatively, several studies have found a trend towards having higher ξ_ion values in galaxies with bluer UV continua, a proxy for galaxies with more dominant young stellar populations (Duncan & Conselice 2015; Bouwens et al. 2015b, 2016). However, Matthee et al. (2017a) show that this is a product of using the UV slope itself (or galaxy SEDs that predominantly trace the rest-UV) as a proxy for dust attenuation. This is also sensitive to the assumed dust model, as can be seen in fig. 4 of Shivaei et al. (2018) where using different dust laws can produce a marginally positive or a negative correlation between ξ_ion and M_UV. Nevertheless, Matthee et al. (2017a) find that ξ_ion increases with decreasing UV luminosity, increasing specific SFR, and increasing H α EW in their sample of z ≈ 2 HAEs.

In Fig. 4, we compare the equivalent width of H α with the derived ξ_ion value for the emission line-selected samples. In addition, we plot the fit to the relation from Tang et al. (2019) derived from z ≈ 1−2 ‘extreme emission line galaxies’ (EELGs), which are selected purely on the basis of high-EW optical emission lines ([O iii] and/or H α). These systems have gas-phase metallicities <0.3 Z_⊙ and mass doubling times on the order of 100 Myr (Maseda et al. 2013, 2014). While these physical properties are broadly similar to those of literature samples of LAEs at higher redshifts, the offset compared with the high-EW LAEs presented here argues for a difference in ξ_ion at a fixed age. This is due to the fact that the EW of Balmer lines scales inversely with the age of the stellar population in these types of systems (e.g. van der Wel et al. 2011). Here and throughout, we specifically refer to the light-weighted age, which is dominated by the most recent generation of star formation: an older (>1 Gyr) stellar population could contribute to the total mass of the galaxy but would not contribute significantly to the UV or H α luminosities. The highest EW objects in the Tang et al. (2019) sample, with ages <10 Myr, have the highest values of ξ _ion. As these galaxies and the high-EW LAEs presented here have little to no dust attenuation (see Section 3.4), any correlation between ξ_ion and age is not driven by using the rest-UV to correct for dust attenuation. At a fixed H α EW, our sample presents a larger value of ξ_ion than the literature LAE/HAE results, implying a lower gas-phase metallicity. However, in all samples a trend with higher ξ_ion at younger ages is found.

Figure 4.

Rest frame H α EW versus ξ_ion for this sample, compared to the literature observations of other line emitters (LAEs and HAEs). The error bars are not independent as a higher H α EW necessarily implies a higher ξ_ion. The relationship between the EW of H α and ξ_ion is consistent between the literature samples, with a strong trend to higher ξ_ion at higher H α EWs. As the H α EW is inversely proportional to the age, this smooth trend suggests that younger ages are the main driver of elevated ξ_ion. Our elevated value at fixed H α EW suggests that these LAEs have a lower (gas phase) metallicity than other samples presented in the literature (see Fig. 5).

Qualitatively, for a constant star formation history, ξ_ion is roughly constant for the first 1 Myr before monotonically declining until the population is approximately 1 Gyr old, when it plateaus again (Fig. 5). This is related to differences in the star formation time-scales probed by H α and the UV, and the effect is strongest between 1 and 10 Myr (or longer when binary stellar evolution is included in the models, e.g. Stanway, Eldridge & Becker 2016). For an instantaneous burst of star formation, the value of ξ_ion drops strongly after 1 Myr. At a fixed age (and hence UV luminosity) for both star formation histories, decreasing metallicity results in an increasing ξ_ion or H α luminosity. So, while young ages result in higher ξ_ion values for the (bursty) HAEs compared to the general population of star-forming galaxies, the observed offset in ξ_ion at fixed H α EW (and hence fixed age) for the high-EW LAEs shown in Fig. 4 suggests that the metallicities must be lower and/or |$f_{\mathrm{esc}}^{\mathrm{ion}}$| must be higher. Independent constraints on the metallicity from, e.g. rest-optical spectroscopy will be required to determine the contribution of each.

Figure 5.

Predicted ξ_ion (left), intrinsic Lyman-α EW (i.e. corrected for the Lyman-α escape fraction; centre), and UV continuum slope (β; right) from the Raiter et al. (2010) stellar population models with constant star formation (solid lines) or an instantaneous burst (dotted lines) at different metallicities (colours), compared to the results presented here. The shaded regions denote our 68 per cent confidence interval in ξ_ion (left), the intrinsic Lyman-α EW for the sample (centre; including the uncertainty on the Lyman-α escape fraction), and β (right), all from the bootstrap iterations described in Section 3.1. We quantify the goodness-of-fit for each combination of star formation history, age, and metallicity in Fig. 6. A colour version is available online.

In Fig. 5, we show the predictions from the Raiter et al. (2010) stellar population models for the evolution of the H α, Lyman-α, and UV luminosities at different metallicities for a constant SFR and an instantaneous burst (their ‘cs5’ and ‘is5’ models, respectively). These models use the stellar tracks, atmospheres, and prescriptions for nebular (line and continuum) emission from Schaerer (2003). For the nebular emission, they assume ionization-bounded nebulae with constant electron temperature and densities, and that all photons are absorbed inside the H ii regions (i.e. |$f_{\mathrm{esc}}^{\mathrm{ion}}$| = 0). We consider models with a Salpeter (1955) initial mass function, with a high-mass cutoff of 500 M_⊙.

At each age and for a given star formation history, the fixed-metallicity models predict our observable quantities: ξ_ion, the EW of Lyman-α, and the UV continuum slope β. We can compare these predicted values to the posterior distributions derived from our bootstrapping procedure and derive a (relative) likelihood for each model. This grid of likelihoods is shown in Fig. 6. While we do not have independent constraints on the star formation history, in both cases the most likely models have metallicities of 0.001 (and ranging between 4 × 10⁻⁴ and 0.008) and ages less than 3 Myr. Lower metallicity models are also permitted with older ages, up to 300 Myr. Although the specific metallicity threshold is dependent on our choice of stellar population model and IMF (cf. Stanway et al. 2016), our median value of ξ_ion is difficult to reproduce with any current set of models even when the high-mass cutoff of the IMF is set to 500 M_⊙ (T. Nanayakkara in preparation). This could be related to the lack of very hard ionizing photons in these models that struggle to reproduce observations of high-z, low-metallicity galaxies (e.g. Kewley, Nicholls & Sutherland 2019; Nanayakkara et al. 2019).

Figure 6.

Scaled likelihood of each of the Raiter et al. (2010) models to reproduce our observations as a function of age and metallicity for a constant star formation history (left) and an instantaneous burst (right). As each combination of metallicity, age, and star formation history produces an estimate for each of our three measured quantities (ξ_ion, Lyman-α EW, and β; Fig. 5); we can calculate a likelihood for that combination by comparing the predictions to our posterior probability distributions for each quantity. The total likelihood is then the product of the three individual probabilities, scaled to a maximum of 1. As shown by the coloured contours, the most likely models have metallicities of 4 × 10⁻⁴ to 0.008 and ages less than 3 Myr, regardless of the star formation history. We cannot, however, conclusively rule out more metal-poor stellar populations with older ages (>10 Myr for constant star formation or >3 Myr for an instantaneous burst). A colour version is available online.

Trainor et al. (2016) determine a maximum ξ_ion value of 10^25.6 Hz erg⁻¹ at 7 per cent Z_⊙ from their sample of z ≈ 2.5 LAEs, some of which have Lyman-α EWs in excess of 100 Å and have UV magnitudes fainter than −18. Based on their measurement of a high-ionization parameter across their full sample, they conclude that age alone does not drive the elevated ξ_ion that they measure and hence these galaxies should be producing ionizing photons in a steady state. Indeed, even at 100 Myr, the models shown in Fig. 5 with metallicities below 3 per cent Z_⊙ are elevated with respect to the canonical value. Metallicity can indeed play a role in having an elevated ξ_ion, as the extreme Lyman-α EWs require metal-poor (<0.02 Z_⊙) stellar populations with ages of 10 Myr or less (Fig. 5 and Ono et al. 2010; Hashimoto et al. 2017a).

The EW of Lyman-α is observed to scale with |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$| (which is correlated with |$f_{\mathrm{esc}}^{\mathrm{ion}}$|⁠, albeit with large scatter Dijkstra, Gronke & Venkatesan 2016), and for our sample of high-EW LAEs, we would expect to have |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$| values close to 1 (Sobral & Matthee 2019). However, we do not observe such large values of |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$|⁠: as shown in Section 3.6, the mean value is 0.217 assuming Case B recombination (cf. Raiter et al. 2010). This matches the value of 0.3 found in Trainor et al. (2015) despite that sample having a mean Lyman-α EW of only 43 Å. Interestingly, recent observations of an EELG at z ≈ 1.8 by Erb et al. (2019) find |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$|to be 0.097, which is also low for its Lyman-α EW. Jaskot et al. (2019) do not find a strong correlation between H α EW, Lyman-α EW, and |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$| in ‘green pea galaxies’, which have extreme ionization fields that are plausibly similar to the galaxies presented here. They note that the highest H α EWs reflect high intrinsic Lyman-α production and hence a large escape fraction is not required to produce the observed EWs of Lyman-α. In fact, the Sobral & Matthee (2019) model predicts the relationship between |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$|and the EW of Lyman-α to flatten with increasing ξ_ion, such that we could expect a similar value of |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$|to the one we measure. The extreme galaxies we select here, then, do not require high values of |$f_{\mathrm{esc}}^{\mathrm{Ly\alpha }}$| to have the very high EW values measured here (259 Å; median). In addition, Smith et al. (2019) find that the escape of Lyman-α photons lags behind the star formation activity by several tens of Myrs; hence, our selection of the youngest star formation episodes via high-EW Lyman-α could preferentially select periods of relatively low escape.

4.1 Alternatives to star formation

So far we have considered only star formation as the source of ionizing photons in these systems. Overall, with the current data, it is difficult to determine what is the main source of ionizing photons inside these galaxies. As previously mentioned, existing models of normal stellar populations can generally reproduce our observables (i.e. Fig. 5). Similar conclusions are drawn by Trainor et al. (2016) for z ≈ 3 LAEs, where they use strong optical emission line ratios to demonstrate that star formation is the dominant source of ionizing photons. In addition, we have demonstrated in Maseda et al. (2018) that the undetected sample plausibly represents the high-EW extension of the distribution of LAEs at these redshifts, a population that is not dominated by active galactic nuclei (AGN).

Thus, we do not believe that AGNs are the only or even the dominant contributor to these systems. However, we cannot conclusively rule out some contribution, either as the primary source in a subset of objects or a secondary source in all objects. Further spectroscopy is required in both the rest frame-UV and rest frame-optical.

5 OUTLOOK AND CONCLUSIONS

In this work, we have identified a population of 35 z ≈ 4–5 high-EW LAEs using deep MUSE spectroscopic data in the UDF, with a median Lyman-α EW of 259 Å. By combining this with ultra-deep Spitzer/IRAC photometry, we have measured H α emission with an EW of 632 Å (210−1600 Å; 68 per cent) from these systems and used it to calculate the ionizing photon production efficiency, ξ_ion (equation 1). Our primary conclusions are as follows:

• Using Spitzer/IRAC photometry to ≈200-h depth and sophisticated techniques for source deblending, we detect H α emission in stacked data for objects with intrinsic UV magnitudes fainter than −16 (Fig. 2).

• The mean value of ξ_ion for these high-EW LAEs is 10^26.28 Hz erg⁻¹ for an escape fraction of zero, with a 68 per cent confidence interval from 10^25.89 to 10^26.56 Hz erg⁻¹ (Section 3.5). At a UV magnitude of −15.7, this value is a factor of 8.7 in excess of the ‘canonical’ value from literature studies of emission line-selected and continuum-selected samples at similar and higher redshifts and also higher than the most extreme values for z ≈ 7 galaxies from Stark et al. (2017) or local galaxies from Chevallard et al. (2018).

• While the values of ξ_ion from the literature are consistent with a constant value at all M_UV (e.g. Lam et al. 2019), our observation lies significantly above this relation (Fig. 3). This naturally follows from selecting objects with bright Lyman-α luminosities (and hence ionizing photon production rate) compared to the UV continuum, i.e. selecting on the EW of Lyman-α.

• Based on the observed trend between H α EW and ξ_ion, as well as models of the evolution of ξ_ion from stellar population synthesis, we determine that an elevated H α luminosity compared with the UV luminosity is likely a natural consequence of a gas-phase metallicity between 4 × 10⁻⁴ and 0.008 and an age younger than 3 Myr, regardless of the star formation history (Figs 5 and 6).

Typical galaxies selected via the Lyman-break technique (e.g. Bouwens et al. 2017) require detections in the rest-UV and hence are biased against the youngest (<10 Myr) galaxies that necessarily have not produced much UV flux. This is true for LAE selections as well (e.g. Lam et al. 2019), even in samples derived from narrow-band imaging (e.g. Harikane et al. 2018), as the relative depth of the emission line detection threshold with respect to the continuum detection threshold prevents these studies from finding the youngest systems. Our MUSE Lyman-α selection, on the other hand, is closer to an H α selection considering it is also a measurement of (emergent) ionizing photons; by requiring non-detections in the rest frame-UV, we preferentially select the youngest systems. With a selection based on maximal line emission and minimal continuum, it is logical that we have found a population of galaxies with elevated values of ξ_ion compared with older, more massive galaxies.

We expect that observing continuum-faint LAEs with lower line EWs would result in a lower determination of ξ_ion as these systems could have older stellar populations, higher metallicities, or both. Indeed, lower-z observations of ξ_ion show significant scatter at low luminosities and stellar masses (M. Paalvast et al. submitted), which is potentially due to a distribution in the ages, star formation histories, and metallicities in this regime. These galaxies are currently unobservable in emission at high-z, but it is possible to measure the product of f_esc and ξ_ion for systems detected in absorption down to M_UV ≈ −16 (Meyer et al. 2019), which could select a population of galaxies that are distinct from LAEs.

Various hydrodynamical simulations predict that the star formation histories of low-mass galaxies in the early Universe are episodic in nature, with periods of star formation followed by more quiescent phases (e.g. Shen et al. 2014; Muratov et al. 2015). Although the precise duty cycle of these star formation episodes is unknown and could vary with properties such as the galaxy halo mass (van der Wel et al. 2011), a single galaxy could undergo multiple episodes of efficient ionizing photon production without significantly building up its stellar mass (e.g. Domínguez et al. 2015). Although a phase with extreme Lyman-α EW may not occur during the lifetime of every galaxy, depending on their star formation and metal enrichment histories, all galaxies should produce an excess of ionizing photons at young ages for several Myrs. Further work is required to fully characterize the duty cycle of intense star formation episodes at high-z, but the steep faint-end slope of the UV luminosity function and the observed high number density of high-EW LAEs imply that a significant number of such events should be taking place across cosmic time: at z = 4–5, high-EW LAEs as selected here represent 6.5 per cent of the full galaxy population at a UV magnitude of −16, based on the Bouwens et al. (2015a) luminosity functions. Such a large cumulative number of these episodes could have a significant contribution to cosmic reionization, especially as LAEs have higher values of |$f_{\mathrm{esc}}^{\mathrm{ion}}$| compared to continuum-selected samples (e.g. Trainor et al. 2015; Marchi et al. 2018) and have number densities at z > 6 that are high enough that they can plausibly contribute significantly to reionization (Drake et al. 2017). As stated above, understanding the distribution of Lyman-α EWs and the relation between LAEs and the general galaxy population at these redshifts is critical, along with the duty cycle, to fully understand the potential contribution of these faint sources to reionization.

A detailed discussion of the implications for reionization is beyond the scope of this work. Many studies combine the measured UV luminosity density ρ_UV, derived from UV luminosity functions, with ξ_ion and |$f_{\mathrm{esc}}^{\mathrm{ion}}$| in order to infer the ionizing emissivity |$\dot{N}_{\mathrm{ion}}$| (e.g. Bouwens et al. 2015b). As we show in Fig. 3, ξ_ion can vary with the UV luminosity (e.g. Duncan & Conselice 2015) and the dust content or metallicity of the galaxy (Shivaei et al. 2018). Additionally, ρ_UV is not well determined down to M_UV < −15 as the majority of the observational constraints come from strong gravitational lensing where the associated systematic uncertainties are large (Bouwens et al. 2017; Atek et al. 2018). Finally, the relationship between LAEs and the general population of galaxies, particularly at these extremely low luminosities, is unclear (e.g. Mesinger et al. 2015).

While these observations are the first step in understanding the ionizing photon production efficiency in such extreme systems, our method necessarily determines the average properties: individual objects could differ significantly. Similar high-EW LAEs at z < 4 could be studied from the ground in the near-IR with detections of H β in extremely long integrations with current instruments, or with future near-IR capabilities provided by the next generation of Extremely Large Telescopes. With the advent of new space-based facilities such as the James Webb Space Telescope (JWST), however, these measurements can be done out to higher redshifts using the brighter H α emission line. Combined with deep rest-UV imaging and spectroscopy to potentially measure |$f_{\mathrm{esc}}^{\mathrm{ion}}$|⁠, we will obtain a more complete picture of the production and escape of ionizing photons from the abundant low-luminosity population of galaxies in the early Universe.

ACKNOWLEDGEMENTS

We would like to thank the anonymous referee for a thoughtful report and suggestions that have improved this manuscript. We are also grateful to everyone involved in the Spitzer Space Telescope mission and everyone at the Spitzer Science Center: we are truly fortunate to have been able to use data from this facility. J. B. acknowledges support by FCT/MCTES through national funds by this grant UID/FIS/04434/2019 and through the Investigador FCT contract no. IF/01654/2014/CP1215/CT0003. S. C. gratefully acknowledges support from Swiss National Science Foundation grant PP00P2|$\_$|163824. We would also like to thank Mauro Stefanon for his assistance with de-blending the IRAC photometry, Pieter van Dokkum for a number of useful suggestions, and Daniel Schaerer for information regarding the stellar population models. This study is based on observations made with ESO telescopes at the La Silla Paranal Observatory under programs IDs 094.A-2089(B), 095.A-0010(A), 096.A-0045(A), and 096.A-0045(B).

Footnotes

1

A 5-σ continuum detection limit as in Maseda et al. (2018) would have resulted in a median rest frame Lyman-α EW of 162 Å.