© ESO 2019

1. Introduction

The reionization of the universe has received significant attention over the last decade. Fundamental unanswered questions remain about both the basic time scale of cosmic reionization and the sources which drive the process. The most obvious sources to power cosmic reionization are star-forming galaxies (e.g., Bouwens et al. 2015; Robertson et al. 2015) and quasars/active galactic nuclei (AGN, Madau & Haardt 2015). One of the reasons why quantifying their respective contribution to cosmic reionization has been challenging is that we cannot detect ionizing photons from these sources directly.

As such, our best estimates for the ionizing emissivity from galaxies or quasars have been based on the emissivity of these sources in the non-ionizing UV-continuum. It has been conventional to convert the galaxy UV luminosity density to an ionizing emissivity (or the rate of ionizing photons per unit volume that reaches the intergalactic medium) ṅ_ion:

(1)

where ρ_UV is the total UV (1500 Å) luminosity density, ξ_ion the ionizing photon production efficiency per unit UV luminosity, and f_esc, UV and f_esc, LyC the fraction of light that is able to escape the galaxy unabsorbed in the non-ionizing UV and ionizing wavelengths, respectively.

Traditionally, the Lyman-continuum photon production efficiency ξ_ion has been estimated by extrapolating from either the rest-frame UV-continuum slope β (Robertson et al. 2013; Duncan & Conselice 2015; Bouwens 2016) or from synthetic stellar population models (Madau et al. 1999). In Bouwens et al. (2016a; hereafter B16), however, it was shown that ξ_ion can be directly inferred from Spitzer/IRAC-based estimates of the Hα fluxes. As demonstrated by Shim et al. (2011) and Stark et al. (2013), observations with Spitzer/IRAC can be used to infer the Hα fluxes as Hα falls in the 3.6 μm (3.8 <  z <  5.0) and 4.5 μm bands (5.1 <  z <  6.6). Given that the vast majority of ionizing photons produced by stars in a galaxy ionize neutral hydrogen which cascades down to produce Hα, the production rate of ionizing photons can be deduced quite straightforwardly from the Hα line flux using quantum mechanics (Leitherer & Heckman 1995).

The availability of two new data sets make it possible to improve upon the estimates in B16, particularly by allowing us to push fainter. The first data set is the Multi Unit Spectroscopic Explorer (MUSE) GTO data over the HUDF and GOODS-S fields (Bacon et al. 2017; Herenz et al. 2017). From these data, it has been possible to construct very large samples of z∼3–6 galaxies, all with spectroscopic redshifts (Inami et al. 2017; Herenz et al. 2017). The advantage of using a purely spectroscopic sample, rather than the photometric samples used in some earlier studies (Smit et al. 2016; Mármol-Queraltó et al. 2016; Rasappu et al. 2016), is that we can derive the Hα fluxes, while ensuring that the Hα line always falls in the desired IRAC filter.

The second of these data sets is the 3.6 μm + 4.5 μm Spitzer/IRAC GREATS observations (Labbé et al. 2014), which has an integration time of 200 h in the deepest parts, resulting in photometric depths at least 0.4 mag deeper than used in the earlier studies on which B16 is based. By taking advantage of the deeper Spitzer/IRAC data and stacking the large numbers of sources we have from MUSE, we can probe the Hα line flux in galaxies which are more representative of the overall population, and have a greater contribution to the overall photon budget, than from the rarest, brightest ones.

In Sect. 2, the observational data are summarized, along with the sample selection criteria. Section 3 describes our procedure for performing photometry, as well as parameter inference. In Sect. 4, we present new measurements of the Lyman-continuum photon production efficiency ξ_ion based on the GREATS data and MUSE spectroscopic redshift samples. Finally, we discuss our results in Sect. 5 and provide a brief summary in Sect. 6. We assume Ω₀ = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹ throughout this paper. All magnitudes are in the AB system (Oke & Gunn 1983).

2. Observational data and sample selection

In this section, we provide a description of all significant data sets that we utilize for our analysis. Table 1 provides a convenient summary.

2.1. Spectroscopic samples

One of the major contributing factors to our being able to push faint in the present study are the large samples of galaxies with spectroscopic redshifts we have from two large MUSE GTO programs: the MUSE-Deep and MUSE-Wide, spectroscopic surveys¹. MUSE-Deep (Bacon et al. 2017) has 116 h of total exposure time over an area of 3.15 × 3.15 arcmin² covering the Hubble eXtreme Deep Field (XDF; Illingworth et al. 2013). MUSE-Deep consists of two components: a shallower 3 × 3 “mosaic” of nine 1 × 1 arcmin² fields that cover the entire MUSE-Deep region, and a deeper, single 1 × 1 arcmin² field named “UDF10” that overlaps with deep ALMA observations (Walter et al. 2016). We use the MUSE-Deep spectroscopic redshift catalog published by Inami et al. (2017), which sources include optical+NIR HST detection as well as blind searches for emission lines in the spectral cube.

The MUSE-Wide survey (Urrutia et al. 2019) complements the MUSE-Deep by surveying a much larger area to probe rare sources. Here we utilize the first data release of MUSE-Wide which consist of 44 1-h pointings (Urrutia et al. 2019). The data cover an area of ≈44 arcmin² distributed over the CANDELS-DEEP WFC3 region (Koekemoer et al. 2011). All spectroscopic redshifts from both MUSE-Deep and MUSE-Wide greater than z = 2.9 are inferred from the Lyman-alpha line.

To make use of our large sample of spectroscopic redshifts to constrain the flux in the Hα line and other lines, we need to segregate sources in redshift windows where a consistent set of strong rest-frame optical emission lines contribute to the IRAC fluxes. The most obvious choices of redshift windows are those where the Hα line lies in either the 3.6 μm or the 4.5 μm band while no other strong lines are present in the other. This criteria is met for sources at 4.532 ≤ z ≤ 4.955 (denoted as z₄) and at 5.103 ≤ z ≤ 5.329 (z₅), as shown in Fig. 1. We define an emission line to be outside of a filter if its central wavelength lies at half of the “minimum in-band transmission” and inside if it lies at the “minimum in-band transmission”. The “minimum in-band transmission” of the IRAC 3.6 μm and 4.5 μm filters are 56.3% and 54.0% respectively. Beyond z = 5.329, the strong emission line [O III] 5006.84 Å enters the 3.6 μm band, while Hα is still present in the 4.5 μm band. Given the potentially large uncertainties that may result in correcting the 3.6 μm flux for the [O III] contribution, we do not consider the 141 sources beyond z = 5.329.

Fig. 1. Illustration of where [S III] 9530.9 Å, Hα, and [O III] 5006.84 Å lie within the Spitzer/IRAC 3.6 μm and 4.5 μm bands at the low and high-redshift ends (left and right panels, respectively) of the redshift intervals selected for this analysis. The transmission curves of the 3.6 μm and 4.5 μm bands are shown as the blue and red lines, respectively. The orange color-shaded regions denote the wavelengths over which relevant lines (Hα and [S III] 9068.6 Å) lie within our defined redshift windows. z1 does not extend beyond z = 2.9 because MUSE sources at z <  2.9 are identified from spectral features other than Lyman α.

The numbers of sources in redshift intervals z₄ and z₅ are small due to their narrow widths. This can be mitigated by extending the same methodology to lower redshifts by taking in to account a number of weaker emission lines. At 3.856 ≤ z ≤ 4.514 (z₃), Hα is present in the 3.6 μm band while the 4.5μm band is contaminated by another potentially strong line [S III] 9068.6 Å. Its strength can be estimated from sources in the 2.900 ≤ z ≤ 3.309 (z₁) and 3.416 ≤ z ≤ 3.829 (z₂) redshift intervals in a manner similar to the measurement of Hα from z₄ and z₅. Although [S III] is present in the 3.6 μm band down to z ≈ 2.5, we do not include sources at redshifts lower than z = 2.9. The reason is that, unlike most MUSE redshifts at z ≥ 2.9, which are measured from the Lyman-alpha line, MUSE redshifts at z <  2.9 are measured from other spectroscopic features (e.g. O II, C III] and absorption lines), which would make our selection more heterogeneous. All redshift ranges and the relevant emission lines are also tabulated in Table 2.

The aforementioned redshift windows are defined only by the major emission lines, i.e. Hα and [O III] 5006.84 Å. The [S III] 9068.6 Å line is also considered to be a major emission line, partially motivated by our observations. Weaker lines including [S III] 9530.9 Å, [N II] 6548.05 Å and 6583.5 Å, and [S II] 6716.0 Å and 6730.0 Å, are considered in our analysis, but are treated as having an amplitude proportional to the stronger lines (see Sect. 4.2).

Figure 2 shows the number of sources with MUSE spectroscopic redshifts that are within our defined redshift windows in MUSE-Deep and MUSE-Wide, respectively. From the MUSE-Deep catalog, we obtained 155, 192, 156, 116, and 36 sources in the redshift intervals z₁, z₂, z₃, z₄, and z₅, respectively. In MUSE-Wide, 131, 146, 164, 84, and 32 sources are identified in these same respective intervals. Note that the number of usable sources becomes lower when we enforce additional selection criteria, such as whether the spectroscopic redshift is consistent with the photometric redshift, and whether accurate photometry is available.

Fig. 2. Number of galaxies vs. spectroscopic redshift based on our MUSE Deep (left) and MUSE wide (right) data sets. The sources highlighted in green (z1) and blue (z2) are those where the [S III] 9068.6 Å line falls in the Spitzer/IRAC 3.6 μm and in the 4.5 μm bands, respectively, and no other strong lines seem likely to be present in the other band. Sources highlighted in cyan (z3) are those where the Hα line is present in the 3.6 μm band and [S III] is in the 4.5 μm band. Sources highlighted in red (z4) and yellow (z5) are those where Hα falls in the 3.6 μm and the 4.5 μm bands, respectively, and no other strong lines are located in the other Spitzer/IRAC band. See Table 2 for the numbers of sources and Fig. 1 for the line location with respect to the IRAC filters.

2.2. IRAC observations

An important aspect to this study is taking advantage of significantly deeper Spitzer/IRAC observations than were available in previous work by Smit et al. (2016) and Rasappu et al. (2016), as discussed in Sect. 1.

We use data from the GOODS Re-ionization Era wide-Area Treasury from Spitzer (GREATS, PI: Labbé) combined with all the previous Spitzer/IRAC observations taken over the CDF-South (Stefanon et al., in prep.). For the 3.6 μm band, the depth within the area we have MUSE spectroscopy ranges from 76 to 278 h, with an average depth of 167 h. For the 4.5 μm band, the depth ranges from 63 to 264 h, with an average depth of 139 h. The full-width-half-max (FWHM) of the point-spread-function (PSF) of the 3.6 μm and 4.5 μm bands are 1.95″ and 2.02″ respectively. By randomly placing 1000 apertures with 0.9″ radius on the background in the area where the MUSE sources are located, and measuring the fluxes within the apertures, we find the 5σ limiting AB magnitudes in the 3.6 μm band and the 4.5 μm band to be 26.6 and 26.5 respectively.

2.3. HST observations

HST observations over our fields allow us to model the SED and provide us with measurements of the UV-continuum flux. For our HST observations over the fields where we have MUSE spectroscopy, we make use of the Hubble Legacy Field (HLF, Illingworth et al. 2016) reduction. Combining observations from 31 observation programs, this HLF reduction constitutes the deepest reduction of the optical and near-infrared data over the GOODS-S region. We use the version 1.5 data² which include the nine filters F435W, F606W, F775W, F814W, F850LP, F105W, F125W, F140W, and F160W. F098M is not utilized due to its insignificant depth over the area where we have MUSE redshift coverage. Limiting AB magnitudes are measured in the same manner as we did for IRAC data, except using smaller apertures of 0.35″ radius. For the typical region within the HLF, the 5σ limiting AB magnitudes we found for F435W, F606W, F775W, F814W, F850LP, F105W, F125W, F140W, and F160W are 28.08, 27.54, 28.11, 26.73, 28.28, 26.59, 26.51, 25.46, and 26.61, respectively. The images we use for our analysis all have a pixel scale of 60 ms.

3. Photometry and parameter inferences

3.1. HST photometry

Our procedure for performing HST photometry can be briefly described as follows. We use the original, non-PSF-matched F850LP images as our detection images. Skipping the PSF-matching step keeps the signal-to-noise at a maximum. Unlike most other studies on high-z galaxy observations, we do not construct our detection image from multiple filters. A combined detection image has boosted S/N but also has slightly different PSF shapes between stars and galaxies due to their color difference. As we will show in Sect. 3.2, the HST detection image has to be used also as the “prior image” for IRAC photometry. Since our sample is limited by detectability in the shallower IRAC data, the benefits of a stable PSF shape are arguably of similar or greater importance to that of an increased S/N in the detection image. We then measure the flux ratios of sources using SExtractor in dual-image mode with the aforementioned detection image and each science image that are PSF-matched using the procedure we describe in Appendix A. We experimented with various combinations of SExtractor parameters and arrived at a set that is optimized for detecting faint and small sources. The optimal parameters are listed in Table B.1. We use the fluxes measured in isophotal apertures as our best estimate of the flux ratios as these apertures give us measurements whose S/N is typically as high or higher than other aperture choices, such as circular and Kron (1980). We then correct all isophotal fluxes for light that falls beyond the isophotal apertures by multiplying with the factor AUTO_F850LP/ISO_F850LP, where AUTO_F850LP is the F850LP flux measured by the Kron aperture. We correct all uncertainties for noise correlation using Eq. (A.20) in Casertano et al. (2000),

(2)

where F_A is the factor which divides the aperture-corrected photometry to give the noise correlation-corrected photometry, s the pixel size, and p the drizzle pixfrac parameter.

We cross-match MUSE sources with HST sources by requiring the presence of a unique HST counterpart within distances of 0.2″ and 0.4″ for MUSE-Deep and MUSE-Wide respectively. We tried cross-matching with maximum allowed offsets ranging from 0.1″ to 0.8″ in steps of 0.1″, and found that the aforementioned values give the highest unique cross-matching rate. When the maximum allowed offset is smaller than 0.2″–0.4″, we lose sources by scatter in astrometry. When the maximum allowed offset is larger than 0.2″–0.4″, we find sources in our MUSE catalogues increasingly matching with more than one source in our HST catalogues, implying that at such large matching radii the identifications are more ambiguous. Our cross-matching rate in MUSE-Deep of 85% is comparable with the fraction of MUSE-Deep sources that have HST counterparts in the UVUDF catalog (Rafelski et al. 2015), which is 89% (Inami et al. 2017).

For MUSE-Wide, our overall cross-matching rate of 76% is slightly lower those cross-matching with the catalogs of Guo et al. (2013) and Skelton et al. (2014), which are 80% and 85% respectively (Urrutia et al. 2019). The deficit is likely caused by the fact that our HST detection scheme is set up for finding high-redshift, faint objects. If we only consider MUSE-Wide objects that have spectroscopic redshifts greater than z = 2.9 (the lowest redshift relevant to this work), then we get a much more consistent cross-matching rate of 55% compared to the 44% and 57% cross-matching rate with Guo et al. (2013) and Skelton et al. (2014) respectively (Urrutia et al. 2019).

3.2. Spitzer/IRAC photometry

Due to the broad PSF of Spitzer/IRAC (FWHM ≈ 2″) and the depth of the GREATS data, source crowding makes it difficult for basic photometric tools (e.g. SExtractor) to make accurate flux measurements for faint sources. Given these challenges, we use the code “MOPHONGO” (Labbé et al. 2015) to model and subtract neighboring sources in the IRAC images within a 12″ × 12″ region surrounding the source. The source morphology in the IRAC images is modeled using PSF-matched HST images, with the total flux of each source being a free parameter. Figure 3 shows an example of the subtraction we obtain of neighboring objects near one of the sources in our samples with a MUSE redshift.

Fig. 3. Illustration of how the subtraction of flux from neighboring sources is performed on the Spitzer/IRAC images to derive accurate photometry for one of our sources in our sample. The “HST data” panel (upper left) shows a cut-out of the original F850LP image centered on a galaxy at z = 4.54. The “HST segmentation” panel (upper right) shows the SExtractor segmentation map over the same region. The “IRAC data” panel (middle left) shows the original Spitzer/IRAC 3.6 μm image. The “model” panel (middle right) shows the best-fit model constructed from segments of psf-matched HST F850LP data. The “residual” panel (lower left) is the residual subtracting the full model from the data. In the “clean” panel (lower right), all objects are subtracted except the target itself. The Spitzer/IRAC 3.6 μm flux is then measured from the “clean” image. Each tile size is 18″ × 18″.

3.3. Absolute UV magnitudes

The absolute UV magnitude, M_UV, is a common metric to quantify the intrinsic luminosity of a source. It allows us to bin objects with similar intrinsic luminosities to obtain meaningful stacked photometry. Following the convention of B16, we measure M_UV at a rest-frame wavelength of 1600 Å. They are computed from the apparent magnitudes in F606W, F775W, F850LP, and F105W for sources in 2.900 <  z <  3.429, 3.429 <  z <  4.250, 4.250 <  z <  5.063, and 5.063 <  z <  5.329 respectively. We correct these magnitudes for the shapes of the transmission curves, and in some cases, the presence of Lyman break (at rest-frame 1216 Å) inside the relevant filter by assuming a simple f_λ ∝ λ⁻² continuum with a sharp Lyman-break cutoff. The left column of Fig. 4 shows the distribution of M_UV at different redshifts.

Fig. 4. Number of sources in each of the spectroscopic redshift intervals as a function of their absolute magnitudes (left column), stellar masses (middle column), and UV continuum slopes (right column). In estimating the strength of the Hα and [S III] lines, we make use of stacks of the Spitzer/IRAC images of the sources from each bin.

3.4. Stellar mass inferences

It is interesting to examine how the emission line strengths and Lyman-continuum photon production efficiency, ξ_ion, depend on the stellar mass. As such, we estimate stellar masses for individual sources to segregate sources into different stellar mass bins. The stellar masses are estimated by fitting synthetic stellar population spectra to the HST+Spitzer photometry with FAST (Kriek et al. 2009). We adopt the stellar population synthesis libraries of Conroy & Gunn (2011), a Chabrier (2003) initial mass function, a delayed exponentially declining star formation history, and a Kriek & Conroy (2013) dust law. Metallicity is allowed to take values of 0.019 (Z_⊙), 0.0096 (≈Z_⊙/2), and 0.0031 (≈Z_⊙/6). Redshift is left as a free parameter because the “photometric redshifts” derived by FAST can be used to strengthen the less-confident spectroscopic redshifts (more on this in Sect. 3.6). Figure 4 shows the distributions of derived stellar masses.

At present, FAST is only capable of modeling the continuum radiation without including contributions from the emission lines. As a result, FAST will try to reproduce the flux in Hα-boosted bands using features in the continuum. This could result in unphysical estimates of the stellar mass. It is useful, therefore, in deriving reliable stellar masses to exclude the Spitzer bands which are contaminated by the Hα emission line, i.e. the 3.6 μm band for sources in the redshift range 3.856 <  z <  4.955, and the 4.5 μm band for sources in the range 5.103 <  z <  5.329.

3.5. UV-continuum slope inferences

We will also want to examine whether the emission line strengths show a dependence on the UV-continuum slope of individual sources – given that the UV-continuum slope β is frequently related to the average light-weighted age of the stellar populations within a galaxy. The UV-continuum slope is parameterized as f_λ ∝ λ^β. Following the convention of Bouwens et al. (2012) and Castellano et al. (2012), we measure β by modeling the photometry of HST filters that correspond to rest-frame wavelengths from 1650 to 2300 Å as a power law. Table C.1 lists the filters used for deriving β for sources at different redshifts. Figure 4 illustrates the distribution of UV continuum slopes derived for our samples.

3.6. Photometric redshift inferences

As a check on the spectroscopic redshifts inferred from the MUSE GTO data, we also derive photometric redshifts for all the sources. The value of doing these checks is clear for the less-certain MUSE redshifts, which are derived from single, low-S/N emission lines, where an asymmetric Ly-α-like profile is not evident. However, the same check is needed even for the more-secure MUSE redshifts (those with high-S/N asymmetric profile of Lyα and/or multiple spectroscopic features). The reason is that some of these secure MUSE sources overlap with brighter, foreground objects which could be confused using our cross-matching procedure. Clearly, such cases have to be excluded from our analysis because the spectroscopic redshifts and the photometry correspond to different objects. We check the spectroscopic redshifts with the photometric redshifts derived with either EAZY (Brammer et al. 2008) or with FAST. We define consistency between spectroscopic redshifts z_spec and photometric redshifts z_phot to be such that |Δ(z_phot − z_spec)| <  0.5.

We derive separate photometric redshifts for each source with FAST and EAZY and only require that one of the two redshift estimates is consistent with a Δz <  0.5 for maximum inclusivity. Fortunately, the two codes give similar results for most sources, despite modest differences in the approaches they utilize.

EAZY derives photometric redshifts by fitting spectral templates, with the flux in one of the bands acting as prior. Although the templates used in EAZY contain emission lines, the line strengths are not free parameters. Since this may have a significant impact on the ability to fit fluxes in bands with Hα line contamination, we excluded those bands from the fits.

After excluding sources with unsatisfactory neighbor-subtraction in the IRAC images (see Sect. 4.1 for the exact criteria), we found that 77% of the sources with less-confident MUSE redshifts are consistent with our photometric redshift estimates. For the 219 sources with confident MUSE redshifts, only two are “incompatible” with our photometric redshift estimates. In both cases, this appears to be due to overlap with foreground sources.

4. Empirical estimate of ξ_ion

4.1. Stacking the IRAC and HST images

For galaxies in each of our subsamples (by UV luminosity, stellar mass, and β), we derive a weighted average flux in each IRAC band. The weight scheme is designed to both maximize the S/N and determine the representative SED. To this end, we stack the IRAC images, after all the neighboring objects have been subtracted, i.e., the “clean” panel in Fig. 3. To optimize the accuracy of the flux measurements we obtain from the stacks, only images with a satisfactory subtraction of their neighbors are considered for inclusion. Quantitatively, we exclude any sources which, in either the 3.6 μm band or the 4.5 μm band, have a reduced chi-squared larger than 1.4, flux uncertainty larger than 30.0 nJy, or contamination by neighbors higher than 550% of the measured F850LP flux. These criteria are set in order to maximize the usable sample size, and are determined by visual inspection of the subset of galaxies at 5.103 <  z <  5.329.

In addition to these cutoffs, each source contributes to the final mean stack based on a three-component weight. Each weight is composed of a reduced chi-squared component, an uncertainty component, and a contamination fraction component,

(3)

where each component is parametrized as follows:

(4)

We examine three different binning schemes in examining the dependencies of ξ_ion on various physical or observational quantities. First, sources are binned according to their absolute UV magnitudes: −21.5 <  M_UV <  −20.5, −20.5 <  M_UV <  −19.5, −19.5 <  M_UV <  −18.5, −18.5 <  M_UV <  −17.5. Second, the binning is done according to the stellar mass: 10.0 >  log₁₀ M_*/M_⊙ >  9.0, 9.0 >  log₁₀ M_*/M_⊙ >  8.0, and 8.0 >  log₁₀ M_*/M_⊙ >  7.0. Third, the binning is performed according to the UV continuum slope: −1.5 >  β >  −2.0, −2.0 >  β >  −2.5, and −2.5 >  β >  −3.1. The distribution of weights w used for the M_UV, M_*, and β-subsample stacks is presented in Fig. D.1.

Of course, for the deep stacks we create for galaxies in different bins of stellar mass M_* and β, we include sources over a very wide range in UV luminosity. For most of these stacks, the brightest and faintest sources typically differ in luminosity by a factor of ≈10, with the most extreme case occurring in a bin with −2.0 >  β >  −2.5 and 3.105 <  z <  3.309. There the brightest source is 1258× more luminous than the faintest source. In that case, the brightest source receives 35% of the total weight for that bin. For the other 18 sources which make up the bin, the contribution is ≈4%.

In addition, in order to provide a larger number of constraints for our modeling, the first three redshift ranges, z₁, z₂, and z₃, are each divided into two pieces for the purpose of the binned stacks. That is, the redshift interval z₁ is divided into two separate intervals 2.900 <  z <  3.105 and 3.105 <  z <  3.309, redshift interval z₂ is divided into 3.416 <  z <  3.623 and 3.623 <  z <  3.829, and redshift interval z₃ is divided into 3.856 <  z <  4.185 and 4.185 <  z <  4.514. Figures 5–7 show the stacks for M_UV, stellar mass, and β respectively.

Fig. 5. Stacked Spitzer/IRAC [3.6] (left) and [4.5] (right) images of sources found in the redshift intervals 2.900 <  z <  3.105, 3.105 <  z <  3.309, 3.416 <  z <  3.623, 3.623 <  z <  3.829, 3.856 <  z <  4.185, 4.185 <  z < 4.514, 4.532 <  z <  4.955, and 5.105 <  z < 5.329. Only sources with satisfactory neighbor subtractions and spectroscopic redshifts consistent with their photometric redshifts are included in the stacks. The stacks are binned into three groups of intrinsic luminosity, −20.5 <  MUV <  −19.5, −19.5 <  MUV <   −18.5, and −18.5 <  MUV <  −17.5. Some neighboring objects remain in the corners because the subtraction is performed within a diameter of 12″ instead of a square of the same size.

Fig. 6. Stacked Spitzer/IRAC [3.6] (left) and [4.5] (right) images of sources vs. their estimated stellar mass in the same redshift intervals considered in Fig. 5. Three different bins in stellar mass are considered: 10.0 <  log10(M*/M⊙)< 9.0, 9.0 <  log10(M*/M⊙)< 8.0, and 8.0 <  log10(M*/M⊙)< 7.0.

Fig. 7. Stacked Spitzer/IRAC [3.6] (left) and [4.5] (right) images of sources vs. their UV-continuum slope β in the same redshift intervals considered in Fig. 5. Three bins in β are considered: −1.5 >  β >  −2.0, −2.0 >  β >  −2.5, and −2.5 >  β >  −3.1.

For symmetry, we have applied an identical weighting scheme when considering quantities derived from the HST photometry, i.e., L_UV and β, for individual sources as we used in weighting these sources in our IRAC stacks. Since our weighting scheme favors sources where the subtraction of neighboring sources is accurate, IRAC-bright objects would not be expected to dominate the stacks.

4.2. Measurement of Hα equivalent widths

To derive the equivalent widths of Hα, we first measure the fluxes on the stacked “clean” IRAC images using circular apertures with a radius of 0.9″. The uncertainties are derived by bootstrapping, which accounts for both source-to-source variations and noise correlation. Specifically, we randomly draw, with replacement, n times from the bin sample, where n is the size of the bin sample. These n sources are then stacked, and the flux measured. We repeat this process 1000 times, and the photometric uncertainty is calculated by taking the standard deviation of the 1000 flux values.

We then fit a simple model spectrum, convolved with the filter transmission curves, to the measured [3.6]−[4.5] colors as a function of redshift. As an example, Fig. 8 shows the measured [3.6 μm]−[4.5 μm] colors for sources in the −20.5 <  M_UV <  −19.5 bin. The model spectrum consists of a power-law continuum (f_λ ∝ λ^β_opt), an Hα emission line, an [S III] 9068.6 Å emission line, and five other secondary emission lines ([N II] 6548.05 Å and 6583.5 Å, [S II] 6716.0 Å and 6730.0 Å, and [S III] 9530.9 Å) whose strengths are fixed relative to that of Hα. Assuming a metallicity of Z = 0.004 = Z_⊙/5, we take the strengths of the secondary lines, relative to Hβ, and the Hα/Hβ ratio from Anders & Fritze-v. Alvensleben (2003) and Leitherer & Heckman (1995), respectively. The non-hydrogen line ratios of Anders & Fritze-v. Alvensleben (2003) are based on the modeling results of Stasińska (1984), who simulate nebular emission given a wide range of physical conditions. Leitherer & Heckman (1995) predict the Hα/Hβ ratio by assuming a gas temperature of 10 000 K, a 10% helium abundance relative to hydrogen, and case B recombination. For simplicity and to avoid overfitting our data, we took the optical continuum slope, the Hα EW, and [S III] EW to be independent of redshift. The fluxes of the secondary lines range from 2% to 7% of the Hα flux. As an example of how our stacked spectra would look, Fig. E.1 shows the implied spectra for 4.532 <  z <  4.955 galaxies in three different bins of UV luminosity M_UV, i.e., −20.5 <  M_UV <  −19.5, −19.5 <  M_UV <  −18.5, and −18.5 <  M_UV <  −17.5.

Fig. 8. [3.6]−[4.5] color as a function of the redshift bin considered for our stacked Spitzer/IRAC images of sources with −20.5 <  MUV <  −19.5. The best-fit color model is shown in blue. The color model is the convolution of the 3.6 μm and 4.5 μm transmission curves with a simple model spectrum, which consists of a power-law continuum, an Hα line, an [S III] 9068.6 Å line, and five secondary lines whose strengths are fixed relative to that of the Hα line (Sect. 4.2).

We find the rest-frame equivalent widths to be EW_Hα = 86 ± 15 Å, EW_Hα = 119 ± 52 Å, and EW_Hα = 327 ± 183 Å for sources in the bins −20.5 <  M_UV <  −19.5, −19.5 <  M_UV <  −18.5, and −18.5 <  M_UV <  −17.5, respectively. Optimal values of all free parameters are listed in Table 3. Similarly, we infer the Hα equivalent widths for M_*-binned and β-binned sources, which are also listed in Table 3.

4.3. Procedure to derive ξ_ion, 0

ξ_ion in Eq. (1) refers to the Lyman-continuum photon production efficiency in the presence of a non-zero f_esc, LyC. Since the precise value of f_esc, LyC is still uncertain, it has become customary (following B16) to leave out this complication by assuming f_esc, LyC = 0,

(5)

The intrinsic production rate of Lyman-continuum photons, which we define as Ṅ(H⁰), is related to the rate of Lyman-continuum photons reaching the intergalactic medium, Ṅ_ion (capital alphabets denote a change from rate densities to rates), by

(6)

where the 1/(1 − f_esc, LyC) factor reflects the fact that Ṅ(H⁰) is derived from the observed, unattenuated Hα flux, which is in turn produced in recombination cascades after unescaped Lyman-continuum photons are absorbed by the neutral hydrogen gas in galaxies.

Using quantum mechanical simulations, Leitherer & Heckman (1995) found the following relation between the Hα luminosity and the intrinsic Lyman-continuum photons production rate:

(7)

We calculate the Hα flux, L(Hα), using the equivalent widths and rest-frame optical slopes derived in Sect. 4.2. The continuum in the Hα-boosted band is calculated by extrapolating from the photometry in the line-free band using the best-fit rest-frame optical slope. The Hα flux is corrected for the expected dust extinction assuming an SMC extinction law, which is (f_esc, UV)^2.6/13.2. Our choice for this particular dust law is motivated by the recent ALMA finding that z ≈ 5−6 galaxies show resemblance to the SMC dust law (Capak et al. 2015; Bouwens et al. 2016b), particularly at lower luminosities (Bouwens et al. 2016b).

Substituting Eqs. (5) and (6) into (1), we can express ξ_ion, 0 as the ratio of the Lyman-continuum photons production rate to the unattenuated rest-frame UV luminosity,

(8)

The observed rest-frame UV luminosity, L_UV, is calculated from M_UV corrected by an escape fraction for UV photons, f_esc, UV, defined by an SMC dust law (Prevot et al. 1984),

(9)

If instead we adopt the Calzetti et al. (2000) dust law, we recover ξ_ion, 0 values that are lower by ≈0.04 dex.

4.4. ξ_ion, 0 vs. M_UV, M_*, and β

For each M_UV, M_*, and β bin, we obtain an optimized value of the Hα EW by fitting the [3.6]−[4.5] colors with respect to redshift (e.g., Fig. 8). While by construction, the L_UV’s of galaxies are fairly uniform across redshift in M_UV bins, larger variations exist in the M_* and β bins. Therefore, we calculate ξ_ion, 0 for each relevant redshift (only z₃, z₄, and z₅, where Hα is imaged in either IRAC bands) and physical property bin, using the same Hα EW derived for that physical property bin, and a distinct L_UV for that redshift bin. The values and uncertainties of ξ_ion, 0 averaged over redshift are listed in Table 3. These average values are calculated by weighting each ξ_ion, 0 value by the inverse of their uncertainty squared.

Figure 9 plots the derived ξ_ion, 0 against M_UV. The inferred values of log₁₀(ξ_ion, 0/Hz erg⁻¹) are , , and Hz erg⁻¹ for −20.5 <  M_UV <  −19.5, −19.5 <  M_UV <  −18.5, and −18.5 <  M_UV <  −17.5, respectively. Given that our selections exclusively include sources with a detectable Lyα emission line, it is possible that our selections to show a higher ξ_ion than the typical star-forming galaxy at z ∼ 3–6 (see Sect. 5.1). We can set a firm lower limit on the population-averaged ξ_ion, 0 by calculating its value in each M_UV bin and then assuming galaxies that are detected in HST images but not by MUSE do not contribute at all to the ionizing photon budget (clearly an extreme assumption). By cross-matching the combined MUSE-Deep and MUSE-Wide catalog with our HST catalog, we found that 11%, 16%, and 16% of the HST-detected sources are detected by MUSE in −20.5 <  M_UV <  −19.5, −19.5 <  M_UV <  −18.5, and −18.5 <  M_UV <  −17.5, respectively. These Lyα fractions imply firm lower limits on population-averaged log₁₀ ξ_ion, 0’s of 24.32, 24.35, and 24.53 for the aforementioned M_UV ranges respectively.

Fig. 9. Present estimates of Lyman continuum photon production efficiency ξion, 0 vs. absolute UV luminosity MUV. The red dashed line denotes the lower limit for the population-averaged ξion, 0 obtained by making the extreme assumption that all sources in UV selections not appearing in our MUSE selections have ξion, 0 = 0 (see Sect. 4.4 for details). For comparison, also shown here are several previous estimates of ξion, 0 by Bouwens et al. (2016a) at z = 3.8–5.0 (light magenta) and z = 5.1–5.4 (dark magenta; using a revised binning, see Sect. 5.1.), Harikane et al. (2018) at z = 4.9 (orange), Matthee et al. (2017) at z = 2.2 (Meurer et al. 1999β-dust correction, green), and Shivaei et al. (2018) at z = 1.4−2.6 (SMC dust correction, blue). The black dashed line denotes the best-fit MUV dependence of ξion, 0 on MUV using the measurements of Bouwens et al. (2016a) and this work. The best-fit relation is ξion, 0 = 0.020(±0.031)(MUV + 20)+25.372(±0.045).

While the redshift evolution of ξ_ion, 0 is an important topic, we do not consider it in this paper due to our moderate number of sources. The lack of a large sample over a range of redshifts forces us to adopt a constant Hα EW model, and so the derived evolution of ξ_ion, 0 is not very meaningful and only reflects the changes in L_UV.

5. Discussion

5.1. Comparison with previous studies

In this paper, we derived ξ_ion, 0 from the Hα equivalent widths measured in broadband Spitzer/IRAC images. This approach was also taken by two previous studies (Bouwens et al. 2016a; Harikane et al. 2018). Comparing with previous work, we find similar values of ξ_ion, 0 in −20.5 <  M_UV <  −18.5 galaxies, and that in the faintest bin (−18.5 <  M_UV <  −17.5) is perhaps ≈0.2 dex higher, but that increase is not significant. Figure 9 plots these ξ_ion, 0 values with respect to M_UV. The spectroscopic sample of Bouwens et al. (2016a) is from Smit et al. (2016) and Rasappu et al. (2016), which, in turn, derive from Ando et al. (2004), Stark et al. (2013, 2017), Balestra et al. (2010), Vanzella et al. (2005, 2006, 2008, 2009). Most of these sources are Lyman-break galaxies photometrically selected using color criteria for spectroscopic follow-up. Our study takes a similar approach. Our MUSE spectroscopic sample is based upon optical+NIR HST detection and blind searches for emission lines in the spectral cubes. Then we verify them with photometric redshifts. As a result, it is unsurprising to find similar ξ_ion, 0’s for comparable M_UV in previous work.

It is interesting to try quantify the dependence of ξ_ion, 0 on M_UV. To examine this, we assume log₁₀ ξ_ion, 0 varies linearly with M_UV and then fit to the present observations and also the spectroscopic sample of Bouwens et al. (2016a). The best-fit line we derive is ξ_ion, 0 = 0.020(±0.031)(M_UV + 20)+25.372(±0.045). While our best-fit favors a slight increase in ξ_ion, 0 towards fainter M_UV’s, the trend is not statistically significant (< 1σ). For the fits, the binning we consider for the z ≈ 5 sample from Bouwens et al. (2016a) has been to include at least 9 sources in each bin, with values of at M_UV = −21.35 and at M_UV = −20.05.

At low redshifts, Shivaei et al. (2018) measured ξ_ion, 0 directly from the Hα line fluxes of 676 galaxies at 1.4 <  z <  2.6 in the MOSDEF spectroscopic survey. By measuring the line flux of an additional Balmer line, Hβ, they are able to accurately correct for the dust attenuation of Hα. They measured a fairly constant log₁₀ ξ_ion, 0 of 25.06 (25.34)³ spanning absolute UV magnitudes from M_UV ≈ −22.5 to M_UV ≈ −19.5. Matthee et al. (2017) calculated ξ_ion, 0 for 588 Hα emitters and 160 Lyα emitters at z = 2.2. Their sample was based on an Hα/Lyα-selection using narrow-band images. They found an increasing trend in ξ_ion, 0 with absolute UV magnitude from M_UV = −23 to M_UV = −17.

5.2. Tentative detection of [S III] 9068.6 Å in the SEDs of z ≈ 2.9–4.5 galaxies

One interesting new finding in our analysis is the detection of [S III] 9068.6 Å in our stack results. It is detected in the stacks of three bins that contain the brightest sources, −20.5 <  M_UV, AB <  −19.5 , 10 <  log₁₀(M_*/M_⊙) < 9, and −1.5 <  β <  −2.0 bins. Our analysis suggests the presence of [S III] 9068.6 Å with a rest-frame EW of ≈110 Å. Its presence can be seen as a “bump” at z ≈ 3.5 in the [3.6]–[4.5] color against redshift (see Fig. 8).

This line is mostly out of the spectral range of spectroscopic surveys at both low and intermediate redshifts which were conducted in the optical, but we can look at theoretical predictions for the line. Anders & Fritze-v. Alvensleben (2003) predicts the flux in the [S III] 9068.6 Å line to be 33% as strong as Hα for galaxies in the metallicity range 0.4–2.5 Z_⊙ and 18% as strong as Hα for galaxies with metallicities of 0.2 Z_⊙. Interestingly, in our brightest magnitude bin, the lines we infer at ≈9000 Å have a measured flux which is ≈11–16% of the flux in the Hα line, which is comparable with predictions of low metallicities.

5.3. Implications for reionization

In evaluating the capacity of star-forming galaxies to drive cosmic reionization, the total ionizing emissivity is typically calculated by multiplying three separate factors: the unattenuated UV luminosity density ρ_UV/f_esc, UV, the ionizing photon production efficiency ξ_ion, and the Lyman-continuum escape fraction f_esc, LyC. In the present analysis, we were able to place constraints on ξ_ion, 0 for galaxies at a time shortly after reionization has completed. If we assume galaxies in the era of reionization (6 <  z <  9) had similar efficiencies in producing ionizing photons as the ones we analyzed here (z ∼ 4–5), we can set limits on the escape fraction of ionizing photons in those galaxies.

Assuming star-forming galaxies drive the reionization of the universe, Bouwens et al. (2015) have shown that the relative escape fraction f_esc, rel ≡ f_esc, LyC/f_esc, UV, and ξ_ion must satisfy the following relation:

(10)

where M_lim is the assumed UV luminosity cut off and f_corr(M_lim) is a correction factor for ρ_UV(z = 8) integrated to different values of M_lim. Bouwens et al. (2015) found that log₁₀(f_corr(M_lim)) could be approximated as 0.02 + 0.078(M_lim + 13)−0.0088(M_lim + 13)². Note that f_corr(M_lim) is close to unity when M_lim = −13, which is a typical limiting magnitude chosen by many studies (e.g., Robertson et al. 2015; Bouwens et al. 2015). The clumping factor, , is commonly chosen to be 3, a value motivated by simulations (e.g., Pawlik et al. 2009).

As indicated by Eq. (10), we have a collective constraint on the product of the efficiency with which star-forming galaxies produce ionizing photons (ξ_ion) and the implied relative escape fraction (f_esc, rel). The product of these two factors cannot be greater than indicated by Eq. (10) or the cosmic reionization would have been completed sooner than observed, i.e., at z ≈ 6.

The average inferred log₁₀ ξ_ion, 0 (in M_UV bins) is 25.36 ± 0.08. Given that faint galaxies provide the dominant contribution to the overall UV luminosity density, we treat all galaxies as having the same ξ_ion, 0 value as the fainter sources we are studying here. If we take the inferred value of ξ_ion, 0 as typical, and approximate ξ_ion = ξ_ion, 0/(1 − f_esc, LyC)≈ξ_ion, 0 for small f_esc, LyC, we determine that the relative escape fraction cannot be larger than ≈8–20% using Eq. (10). This relation and the constraint on the relative escape fraction are visually presented in Fig. 10.

Fig. 10. Relative escape fraction fesc, rel ≡ fesc, LyC/fesc, UV as a function of ξion, 0 is plotted in green, assuming star-forming galaxies solely drove reionization. Since the dust attenuation estimated for (non-ionizing) UV wavelengths is negligible (fesc, UV = 1), fesc, rel equals fesc, LyC. We also approximated that, for small fesc, LyC, ξion ≈ ξion, 0. The green region was derived in the analysis of Bouwens et al. (2015) by making use of important constraints on reionization (see Table 1 of Bouwens et al. 2015 for a complete list). The blue region shows the range of ξion, 0 inferred in the present study for faint galaxies, i.e. ξion, 0 = 25.36 ± 0.08. The corresponding constraints we can place on fesc, rel of ≈8–20% are indicated by the dashed lines. See Sect. 5.3.

The average UV-continuum slope of sources within −20.5 <  M_UV <  −17.5 at 3.856 <  z <  5.329 is β = −2.26. This means the correction for dust extinction required is negligible (f_esc, UV = 1) according to the SMC dust law (Eq. (9)). In this case, the relative escape fraction, f_esc, rel, equals the Lyman-continuum escape fraction, f_esc, LyC. Our constraints on f_esc, rel (or f_esc, LyC) is consistent with the new Lyman-continuum escape fraction results from Steidel et al. (2018), which imply f_esc, LyC values of ≈0.09 for sub-L* and brighter galaxies.

6. Summary

In this paper, we measured the EWs of Hα and the Lyman-continuum photon production efficiency ξ_ion, 0 for galaxies fainter than 0.2 L^* in the redshift interval z ∼ 3–5. Because Hα is a recombination line, its EW provides a useful measurement of ξ_ion, 0, the intrinsic rate at which ionizing photons are produced per UV-continuum photon. Since faint galaxies likely dominated the ionizing photon budget, measurements of ξ_ion, 0 for faint galaxies allow us to better understand the role of star-forming galaxies played in cosmic reionization.

We are able to extend ξ_ion, 0 measurements to uniquely low luminosities thanks to many spetroscopic measurements for faint sources from the MUSE GTO program and the 200-h Spitzer/IRAC data now available from GREATS + other Spitzer/IRAC programs. This combined data set constitutes the deepest available Spitzer/IRAC imaging over any part of the sky. To measure accurate IRAC photometry, we use the deep HST images as priors for modeling the surface brightness profiles of sources in IRAC images using the code MOPHONGO.

The large number of spectroscopic redshifts we have available from the MUSE-Deep and MUSE-Wide programs allow us to segregate sources into a few distinct redshift intervals where specific strong emission lines fall within distinct Spitzer/IRAC filters (see Fig. 1 and Table 2). From the redshift dependence of the [3.6]−[4.5] colors we can estimate the Hα EW, and thus infer ξ_ion, 0.

Sources are subdivided into different subsamples according to their physical properties (M_UV, M_*, and β), and are stacked (see Eqs. (3) and (4)). To account for the potential contamination of [S III] 9068.6 Å, we consider galaxies at lower redshifts that give us leverage on the [S III] equivalent width in the same manner as z ≈ 4–5 galaxies give on Hα. We fit the observed [3.6]–[4.5] colors across z ≈ 2.9 to z ≈ 5.3 to constrain, simultaneously, the Hα EW, the [S III] EW, and the rest-frame optical continuum slope. We measure a rest-frame Hα EW of 403–2818 Å and a rest-frame [S III] EW of ≈110 Å.

From our inferred Hα EWs, we estimate an average log₁₀ ξ_ion, 0 of 25.36 ± 0.08 for sources between −20.5 <  M_UV <  −17.5. As such, we have been able to estimate ξ_ion, 0 for UV luminosities 1.5 mag fainter than was possible in Bouwens et al. (2016a), and 2.5 mag fainter than Harikane et al. (2018). Combining our new results with those from previous studies that probe brighter magnitudes, we do not find any statistically significant (< 1σ) dependence of ξ_ion, 0 on the UV luminosity of star-forming galaxies. The larger uncertainties at high redshifts, however, do not imply inconsistency with the trend found at lower redshifts (Matthee et al. 2017). If we take our derived ξ_ion, 0’s as typical, they imply a relative escape fraction no higher than ≈8–20% for faint galaxies.

Our ξ_ion, 0 values are potentially biased high due to the selection of sources which show Lyα in emission in the MUSE data, potentially selecting those sources going through an active burst. To control for this potential bias, we will investigate how the Lyman-continuum photon production efficiency depends on the EW of the Lyα emission line in a forthcoming paper and combine with the Lyα escape fraction.

Acknowledgments

We acknowledge useful discussions with Jorryt Matthee and Renske Smit, and also support from NASA grant NAG5-7697, NASA grant HST-GO-11563, and NWO vrij competitie grant 600.065.140.11N211. JR acknowledges support from the ERC starting grant 336736-CALENDS.

References

Appendix A: PSF matching across the various HST bands

We ran SExtractor on the F850LP (“z”)-band image with a set of parameters optimized for detecting stars that have intermediate brightness, as listed in Table B.2. We select stars with CLASS_STAR greater than or equal to 0.1, and AB magnitudes within 16.0 ≤ z_850, AB ≤ 20.0. We choose the size of the PSFs to be 91 × 91 pixels, or 5.46″ × 5.46″. Stars close to the edge (distance from the star’s center to the edge less than 46 pixels) are excluded. We also require a minimum separation 1.365″ (a quarter of the PSF size) between stars and pixels belonging to other sources. There are 43 stars in total within the MUSE-Deep and MUSE-Wide regions that meet these criteria. For each of the HSF bandpasses, PSFs are created by weighting, shifting, normalizing, and median stacking cut-outs of the 43 stars.

Our PSF-matching kernels (from a given HST band to the F160W band) are constructed using the create_matching_kernel Python module from photutils. Noise in the Fourier transforms is filtered with a “split cosine bell window”, which has a parameter controlling the percentage of tapered values, and another controlling the fraction of the array size at which tapering begins. We search for the best kernel by fitting the reproduced radial profile of F160W PSF to the observed one. Figure A.1 shows the radial profiles of reproduced F160W PSFs. It is clear that the PSF-matched kernels are accurate to ≲1% to a radius of 2.5″.

Fig. A.1. Upper sub-panels: comparison of the circularly-averaged radial profile in the F160W band (black curves) with similar radial profiles of the PSFs in the other bands convolved with the corresponding best-fit kernels (blue curves). Lower sub-panels: fractional residual between the profiles is shown. The gray region denotes ±5%.

Appendix B: SExtractor parameters

In Tables B.1 and B.2, we present the SExtractor parameters we use for color measurements and the detection of stars.

Appendix C: Derivation of UV continuum slope

Appendix D: Distribution of the weights utilized in our stack results

In coadding the “cleaned” IRAC images from many sources to create deep stacks, we weight images in the stacks according to the quality of the stacks, contamination from the neighbors, and estimated flux uncertainties in the stack (see Sect. 4.1). To illustrate the range of weights (i.e., w = w_χ² + w_err + w_cont) we use in stacking sources as a function of M_UV, M_*, and β, we show histograms of the weight distributions in Fig. D.1.

Fig. D.1. Distribution of total weights used for individual sources in the MUV-, M*-, and β-binned stacks (see Sect. 4.1).

Appendix E: Illustration of the best-fit spectra implied for different IRAC stacks in our analysis

In deriving the implied Hα EWs for sources, we make use of deep stacks of the Spitzer/IRAC observations from a large number of sources. To provide our readers of how the SED stacks for individual sources look, we have included Fig. E.1 showing the IRAC stack results for galaxies in the redshift range 4.532 <  z <  4.955 for three different bins in UV luminosity −20.5 <  M_UV <  −19.5, −19.5 <  M_UV <  −18.5, and −18.5 <  M_UV <  −17.5.

Fig. E.1. Illustration of our stacked photometric measurements (solid circles) in the 3.6 μm and 4.5 μm bands for galaxies at 4.532 <  z <  4.955 in three different bins of UV luminosity MUV −20.5 <  MUV <  −19.5 (left), −19.5 <  MUV <  −18.5 (center), and −18.5 <  MUV <  −17.5 (right). The width of the 3.6 μm and 4.5 μm bands are shown with horizontal error bars, while the vertical error bars are 1σ. The lines show example model spectra using the line ratios assumed in our analysis (Sect. 4.2)

Appendix F: Consistency check with Smit et al. (2016)

As a test on the consistency of the Hα luminosities inferred in this study with other work, we consider the Smit et al. (2016) study which derived Hα luminosities for many sources over the GOODS North and South fields and conduct a comparison. We compute the Hαluminosity from the color difference:

where ([3.6]–[4.5])_cont = 0.04 is the average continuum color derived from bright (−20.8 <  M_UV <  −19.4) sources at 2.514 <  z <  3.829, and assume f_3.6/f_4.5 equals (f_3.6/f_4.5)_cont. Smit et al. (2016) compute the Hα luminosities from the difference between the observed f_3.6 and the f_3.6 expected from a best-fit SED without emission lines. Table F.1 compares our derived Hα luminosities with those of the overlapping sources in Smit et al. (2016). We find good agreement with Smit et al. (2016) within 1σ for most sources.

All Tables

All Figures

Fig. 1. Illustration of where [S III] 9530.9 Å, Hα, and [O III] 5006.84 Å lie within the Spitzer/IRAC 3.6 μm and 4.5 μm bands at the low and high-redshift ends (left and right panels, respectively) of the redshift intervals selected for this analysis. The transmission curves of the 3.6 μm and 4.5 μm bands are shown as the blue and red lines, respectively. The orange color-shaded regions denote the wavelengths over which relevant lines (Hα and [S III] 9068.6 Å) lie within our defined redshift windows. z1 does not extend beyond z = 2.9 because MUSE sources at z <  2.9 are identified from spectral features other than Lyman α.

In the text

Fig. 2. Number of galaxies vs. spectroscopic redshift based on our MUSE Deep (left) and MUSE wide (right) data sets. The sources highlighted in green (z1) and blue (z2) are those where the [S III] 9068.6 Å line falls in the Spitzer/IRAC 3.6 μm and in the 4.5 μm bands, respectively, and no other strong lines seem likely to be present in the other band. Sources highlighted in cyan (z3) are those where the Hα line is present in the 3.6 μm band and [S III] is in the 4.5 μm band. Sources highlighted in red (z4) and yellow (z5) are those where Hα falls in the 3.6 μm and the 4.5 μm bands, respectively, and no other strong lines are located in the other Spitzer/IRAC band. See Table 2 for the numbers of sources and Fig. 1 for the line location with respect to the IRAC filters.

In the text

Fig. 3. Illustration of how the subtraction of flux from neighboring sources is performed on the Spitzer/IRAC images to derive accurate photometry for one of our sources in our sample. The “HST data” panel (upper left) shows a cut-out of the original F850LP image centered on a galaxy at z = 4.54. The “HST segmentation” panel (upper right) shows the SExtractor segmentation map over the same region. The “IRAC data” panel (middle left) shows the original Spitzer/IRAC 3.6 μm image. The “model” panel (middle right) shows the best-fit model constructed from segments of psf-matched HST F850LP data. The “residual” panel (lower left) is the residual subtracting the full model from the data. In the “clean” panel (lower right), all objects are subtracted except the target itself. The Spitzer/IRAC 3.6 μm flux is then measured from the “clean” image. Each tile size is 18″ × 18″.

In the text

	Fig. 4. Number of sources in each of the spectroscopic redshift intervals as a function of their absolute magnitudes (left column), stellar masses (middle column), and UV continuum slopes (right column). In estimating the strength of the Hα and [S III] lines, we make use of stacks of the Spitzer/IRAC images of the sources from each bin.
In the text

Fig. 5. Stacked Spitzer/IRAC [3.6] (left) and [4.5] (right) images of sources found in the redshift intervals 2.900 <  z <  3.105, 3.105 <  z <  3.309, 3.416 <  z <  3.623, 3.623 <  z <  3.829, 3.856 <  z <  4.185, 4.185 <  z < 4.514, 4.532 <  z <  4.955, and 5.105 <  z < 5.329. Only sources with satisfactory neighbor subtractions and spectroscopic redshifts consistent with their photometric redshifts are included in the stacks. The stacks are binned into three groups of intrinsic luminosity, −20.5 <  MUV <  −19.5, −19.5 <  MUV <   −18.5, and −18.5 <  MUV <  −17.5. Some neighboring objects remain in the corners because the subtraction is performed within a diameter of 12″ instead of a square of the same size.

In the text

	Fig. 6. Stacked Spitzer/IRAC [3.6] (left) and [4.5] (right) images of sources vs. their estimated stellar mass in the same redshift intervals considered in Fig. 5. Three different bins in stellar mass are considered: 10.0 <  log10(M*/M⊙)< 9.0, 9.0 <  log10(M*/M⊙)< 8.0, and 8.0 <  log10(M*/M⊙)< 7.0.
In the text

	Fig. 7. Stacked Spitzer/IRAC [3.6] (left) and [4.5] (right) images of sources vs. their UV-continuum slope β in the same redshift intervals considered in Fig. 5. Three bins in β are considered: −1.5 >  β >  −2.0, −2.0 >  β >  −2.5, and −2.5 >  β >  −3.1.
In the text

Fig. 8. [3.6]−[4.5] color as a function of the redshift bin considered for our stacked Spitzer/IRAC images of sources with −20.5 <  MUV <  −19.5. The best-fit color model is shown in blue. The color model is the convolution of the 3.6 μm and 4.5 μm transmission curves with a simple model spectrum, which consists of a power-law continuum, an Hα line, an [S III] 9068.6 Å line, and five secondary lines whose strengths are fixed relative to that of the Hα line (Sect. 4.2).

In the text

Fig. 9. Present estimates of Lyman continuum photon production efficiency ξion, 0 vs. absolute UV luminosity MUV. The red dashed line denotes the lower limit for the population-averaged ξion, 0 obtained by making the extreme assumption that all sources in UV selections not appearing in our MUSE selections have ξion, 0 = 0 (see Sect. 4.4 for details). For comparison, also shown here are several previous estimates of ξion, 0 by Bouwens et al. (2016a) at z = 3.8–5.0 (light magenta) and z = 5.1–5.4 (dark magenta; using a revised binning, see Sect. 5.1.), Harikane et al. (2018) at z = 4.9 (orange), Matthee et al. (2017) at z = 2.2 (Meurer et al. 1999β-dust correction, green), and Shivaei et al. (2018) at z = 1.4−2.6 (SMC dust correction, blue). The black dashed line denotes the best-fit MUV dependence of ξion, 0 on MUV using the measurements of Bouwens et al. (2016a) and this work. The best-fit relation is ξion, 0 = 0.020(±0.031)(MUV + 20)+25.372(±0.045).

In the text

Fig. 10. Relative escape fraction fesc, rel ≡ fesc, LyC/fesc, UV as a function of ξion, 0 is plotted in green, assuming star-forming galaxies solely drove reionization. Since the dust attenuation estimated for (non-ionizing) UV wavelengths is negligible (fesc, UV = 1), fesc, rel equals fesc, LyC. We also approximated that, for small fesc, LyC, ξion ≈ ξion, 0. The green region was derived in the analysis of Bouwens et al. (2015) by making use of important constraints on reionization (see Table 1 of Bouwens et al. 2015 for a complete list). The blue region shows the range of ξion, 0 inferred in the present study for faint galaxies, i.e. ξion, 0 = 25.36 ± 0.08. The corresponding constraints we can place on fesc, rel of ≈8–20% are indicated by the dashed lines. See Sect. 5.3.

In the text

	Fig. A.1. Upper sub-panels: comparison of the circularly-averaged radial profile in the F160W band (black curves) with similar radial profiles of the PSFs in the other bands convolved with the corresponding best-fit kernels (blue curves). Lower sub-panels: fractional residual between the profiles is shown. The gray region denotes ±5%.
In the text

	Fig. D.1. Distribution of total weights used for individual sources in the MUV-, M*-, and β-binned stacks (see Sect. 4.1).
In the text

Fig. E.1. Illustration of our stacked photometric measurements (solid circles) in the 3.6 μm and 4.5 μm bands for galaxies at 4.532 <  z <  4.955 in three different bins of UV luminosity MUV −20.5 <  MUV <  −19.5 (left), −19.5 <  MUV <  −18.5 (center), and −18.5 <  MUV <  −17.5 (right). The width of the 3.6 μm and 4.5 μm bands are shown with horizontal error bars, while the vertical error bars are 1σ. The lines show example model spectra using the line ratios assumed in our analysis (Sect. 4.2)

In the text