THE SPECTRAL EVOLUTION OF THE FIRST GALAXIES. II. SPECTRAL SIGNATURES OF LYMAN CONTINUUM LEAKAGE FROM GALAXIES IN THE REIONIZATION EPOCH

The SED model results presented in this paper are based on the Yggdrasil spectral synthesis code (Zackrisson et al. 2011b), which mixes single-age stellar population spectra to simulate arbitrary star formation histories. These are used as inputs in the photoionization code Cloudy (Ferland et al. 1998). The result is an age sequence of model SEDs including stellar continuum emission, nebular continuum emission, and nebular emission lines. In this paper, we adopt Starburst99 stellar SEDs (Leitherer et al. 1999) generated using Padova-AGB stellar evolutionary tracks (Vázquez & Leitherer 2005) at metallicities Z = 0.0004–0.020, assuming a Kroupa (2001) universal stellar initial mass function (IMF) over the mass range 0.1–100 M_☉. Since we are focusing on galaxies well into the reionization epoch, we limit the discussion to objects at z ≳ 7, with ages ≲ 7 × 10⁸ yr (roughly the age of the universe at z = 7). Population III galaxies (at zero or near-zero metallicities) may also exist at high redshifts (e.g., Stiavelli & Trenti 2010) and are predicted to exhibit pronounced LyC leakage (Johnson et al. 2009). While such objects likely have too low masses to contribute substantially to the reionization of the universe, we briefly discuss models for such galaxies in Section 4.3 using stellar SEDs from Raiter et al. (2010).

Yggdrasil assumes that the nebular component of the overall galaxy SED can be treated as originating from a single, spherical, and isotropic H ii region (schematically depicted in Figure 1) with constant density and filling factor. In the case of radiation-bounded nebulae with holes (Figure 1(a)), model SEDs with different f_esc are generated by varying the covering factor of the nebula in Cloudy. In the case of density-bounded nebulae (Figure 1(b)), different f_esc are instead achieved by truncating the nebula at a fixed fraction of the theoretical Strömgren radius. In this paper, we adopt a constant hydrogen number density n_H = 100 cm⁻³, a filling factor of 0.01, and a gas metallicity identical to that of the stars. However, the primary diagnostics that we discuss (Balmer recombination lines and nebular continuum in the rest-frame UV/optical) are not sensitive to these parameter choices. We have, moreover, checked that our predictions for radiation-bounded nebulae with holes (Figure 1(a)) are in reasonable agreement with those of the Inoue (2010, 2011) model, which is based on a slightly different computational machinery.

Some of the ionizing photons produced by massive, hot stars are absorbed by the neutral hydrogen in the ISM and produce prominent H ii regions. As a result, photons at rest-frame wavelengths λ < 912 Å get reprocessed into emission lines and nebular continuum flux at longer wavelengths. A substantial fraction of the observed rest-frame UV/optical fluxes of young and/or star-forming galaxies is therefore expected to come from nebular emission rather than direct star light (e.g., Zackrisson et al. 2008; Schaerer & de Barros 2009; Inoue 2011). Since the relative contribution from nebular emission to the overall galaxy spectrum becomes smaller if some fraction of the ionizing photons escape directly into the IGM, information about the LyC escape fraction is imprinted in the non-ionizing (λ > 912 Å) part of the SED. This is demonstrated in Figure 2, where we present model SEDs for young starburst galaxies (age 10 Myr) at metallicities Z = 0.020 (Figure 2(a)) and Z = 0.0004 (Figure 2(b)). The line colors in each panel represent different f_esc (0.0, 0.3, 0.5, and 0.7). At both metallicities, lower f_esc values imply stronger emission lines and a more prominent Balmer jump at 3646 Å. In Figure 2(b), the entire slope of the UV continuum changes with f_esc, because of the higher relative contribution from nebular emission at low metallicity. Similar changes in the UV slope are seen at higher metallicities as well, albeit at younger starburst ages. The noticeable dependence of several spectral features on f_esc suggest that it should be possible to assess the LyC escape fraction from galaxy SEDs without actually measuring the LyC flux. The question is simply how to best retrieve the f_esc information from the rest-frame UV/optical SED that can actually be observed in the reionization epoch (λ > 1216 Å, since the flux at shorter wavelengths is absorbed by the neutral IGM).

Zoom In Zoom Out Reset image size

The power-law slope β (f_λ∝λ^β) of the UV continuum f_λ can in principle be used to gauge f_esc, since the addition of nebular emission to a young, purely stellar SED tends to shift the slope of the UV continuum in the redward direction and increase β (e.g., Raiter et al. 2010; Robertson et al. 2010). The interpretation becomes ambiguous, however, unless a very blue UV slope is detected (β ≲ −3). The distribution of β slopes among galaxies in the reionization epoch is still a matter of debate. Bouwens et al. (2010) reported an average $\overline{\beta }\approx -3$ for z ≈ 7, but later studies have given no support to a typical slope this extreme and instead indicate $\overline{\beta }\approx -2$ (e.g., Schaerer & de Barros 2010; McLure et al. 2010; Bouwens et al. 2012b; Finkelstein et al. 2012b; Dunlop et al. 2013; Rogers et al. 2013). It is possible, however, that a subset of objects still display UV slopes substantially bluer than this (e.g., Finkelstein et al. 2012b; Jiang et al. 2013).

The problem in inferring f_esc from β is that β also depends on age, metallicity, and dust reddening (e.g., Schaerer & Pelló 2005; Bouwens et al. 2010). This is demonstrated in Figure 3, where we plot the age dependence of β for Z = 0.020 (Figure 3(a)) and Z = 0.0004 (Figure 3(b)) galaxies with either single-age stellar populations (dashed lines) or constant star formation rates (SFRs; solid lines). Different line colors are used to indicate the impact of the LyC escape fraction in the case of a radiation-bounded nebula with holes (Figure 1(a)).

Zoom In Zoom Out Reset image size

The arrow shows the shift in β (Δ(β) ≈ 0.2) predicted for A(V) = 0.2 mag of V-band extinction of the stellar continuum (corresponds to ≈0.6 mag at a rest wavelength of 1500 Å) in the case of f_esc = 0, assuming the Cartledge et al. (2005) average Large Magellanic Cloud (LMC) extinction curve. This amount of UV/optical extinction is consistent with several current estimates of the amount of dust affecting the UV slopes of z ≳ 7 galaxies (e.g., Finkelstein et al. 2012b; Dunlop et al. 2013; Bouwens et al. 2012b; Wilkins et al. 2013). At fixed A(V) = 0.2 mag, the Calzetti et al. (2000) attenuation curve gives a similar result (Δ(β) ≈ 0.2). As seen in Figure 3, the effect of f_esc on β is very small for metal-rich galaxies (Figure 3(a)) but becomes substantial at low metallicities (Figure 3(b)) due to the greater importance of nebular emission in the latter case. However, even at metallicities as low as Z = 0.0004, it is only for extremely blue UV slopes (β ⩽ −2.7) that any sort of constraint on f_esc can be set from β alone.

Here, we use the original definition of β given by Calzetti et al. (1994), in which the UV slope is derived from the overall continuum flux (stellar and nebular) in 10 wavelength intervals (chosen to avoid stellar and interstellar absorption features) in the range ≈1270–2580 Å. The exact value of β at any given age depends a bit on the exact wavelength region over which the slope is measured (Raiter et al. 2010), but the overall trends seen in Figure 3 remain the same, regardless of the definition of β used. It should be noted, however, that the 2175 Å dust feature can interfere with a few of the Calzetti et al. (1994) bands and make the UV slope behave in unexpected ways. In this paper, we therefore refrain from discussing attenuation receipts with very prominent 2175 Å bumps, like the average Milky Way extinction curve (e.g., Gordon et al. 2009).

We suggest that by combining the UV slope β with the EW of a Balmer line such as Hβ, it should be possible to identify galaxies with high escape fractions (f_esc ⩾ 0.5) as long as the UV slope is β ⩽ −2.3 after dust reddening corrections (not far from the typical slope measured at z ≳ 6, which is likely to be at least slightly affected by dust; McLure et al. 2010; Bouwens et al. 2012b; Finkelstein et al. 2012b; Dunlop et al. 2013; Wilkins et al. 2013). This idea is demonstrated in Figure 4, where we plot β against the rest-frame EW(Hβ) for Z = 0.020 for constant-SFR, dust-free models with different LyC escape fractions. The results are shown for both geometries depicted in Figure 1, i.e., a radiation-bounded nebula with holes (solid lines) and a density-bounded nebula (dashed lines).

Zoom In Zoom Out Reset image size

The high-leakage models (f_esc ⩾ 0.5; orange, green, and blue lines) are all shifted to the right in Figure 4 (i.e., toward lower EW(Hβ)) compared with the no-leakage case (red lines). Since the model tracks remain separated up to the highest age considered (700 Myr; approximately the age of the universe at z = 7), this diagram suggests that f_esc can be assessed from a simultaneous measurement of EW(Hβ) and β (additional diagnostics can, however, improve the quality of this estimate, as discussed below). The results for radiation-bounded and the density-bounded models are next to identical at f_esc > 0 (and, by definition, completely identical for f_esc = 0). In the following, we will therefore focus on the radiation-bounded models.

At this metallicity (Z = 0.020), β changes very little as f_esc is increased, simply due to the fact that nebular emission has such a small impact on the UV slope (Figure 3). EW(Hβ) of course drops with increasing f_esc, but since nebular emission produces both the Hβ emission line and part of the continuum below the line, the decrease is not as dramatic as one would expect if one assumed the underlying continuum to be purely stellar. Similar diagnostic diagrams can in principle be formed for other Balmer lines (Hα, Hγ, Hδ, etc.), but Hβ seems to be a good compromise between observational limitations (Hα redshifts out of the JWST/NIRSpec range already at z ⩾ 6.6) and expected diagnostic value (Balmer lines beyond Hβ are weaker and have lower EWs). Other definitions of β give rise to slightly different diagnostic diagrams, but the overall trends remain the same. While the [O iii]λ5007 emission line can be measured over approximately the same redshift interval as Hβ, and will in many cases be stronger, the EW of this line is not as suitable as a diagnostic, since it is very sensitive to the ionization parameter of the gas.

Since metallicity is an important parameter in regulating the impact of nebular emission in galaxies (Zackrisson et al. 2008), one expects diagnostic diagrams such as Figure 4 to show some metallicity dependence. In Figure 5, we display the same Z = 0.020, constant-SFR model tracks for radiation-bounded nebulae as in Figure 4, alongside otherwise identical models at Z = 0.004 and Z = 0.0004. As expected, β evolves more strongly with f_esc at low metallicities (see Figure 3), whereas EW(Hβ) displays a milder metallicity dependence. The latter effect is due to the greater contribution from nebular emission to both the Hβ emission line and the underlying continuum at low metallicities. While a rough assessment of f_esc can be obtained even if the metallicity is unknown, a more precise estimate requires some handle on the metallicity. This is most readily obtained from the strengths of emission lines like [Ne iii]λ3869, [O ii]λ3727, and [O iii]λ4959, 5007 (e.g., Nagao et al. 2006), since data on these lines come "for free" when doing spectroscopy over the rest-frame wavelength interval (≈1200–5000 Å) required to measure β and EW(Hβ). It should be noted, however, that many of the metallicity calibrations developed for these lines are based on the assumption of no LyC leakage and may potentially give biased results if blindly applied to density-bounded nebulae with high f_esc (e.g., Giammanco et al. 2005).

Zoom In Zoom Out Reset image size

As illustrated in Figure 6, temporal variations of the SFR can affect both EW(Hβ) and β and need to be considered when attempting to assess the LyC escape fraction. Galaxies with increasing SFRs simply remain longer in the high EW(Hβ), blue β part of the EW(Hβ)–β diagram, which makes the determination of f_esc easier than in the case of a constant SFR (Figure 4). This is shown in Figure 6(a), where Z = 0.020 models with exponentially increasing SFRs (SFR(t)∝exp (− t/τ) with τ = −10⁸ yr) at various f_esc are used. Models with decreasing SFRs are, however, more troublesome. In Figure 6(b), we show the behavior of models with exponentially declining SFRs (τ = 10⁸ yr). These models evolve faster than the constant-SFR models of Figure 4 and reach slightly redder β values at the highest ages considered. However, as there is no overlap between the different f_esc tracks, it should still be possible to estimate the LyC escape fraction. In Figure 6(c), we consider a more extreme model with a very short burst (10⁷ yr) of constant star formation, after which the SFR immediately drops to zero. Since the LyC flux decreases very quickly once star formation has ceased, models of this type evolve sharply to the right in the EW(Hβ)–β diagram once the burst is over. During the post-burst phase, there is a brief period where model lines corresponding to different f_esc overlap, which means that an accurate assessment of f_esc would be very difficult in this part of the diagram. For these particular models, this happens at ages ≈13–15 Myr (β ≈ −2.4), when the LyC flux has dropped to 2%–15% of what it was right before star formation ceased. We stress, however, that galaxies will only rarely be caught in this brief transition phase.

Zoom In Zoom Out Reset image size

Single-component models of the type used in Figures 6(a)–(c) are often assumed to be adequate in the interpretation of observational data of high-redshift galaxies, but real galaxies may well exhibit more complicated star formation histories with several bursts with intermittent periods of low star formation in between. In Figure 6(d), we consider a sequence of equal-strength, 10⁷ yr long, constant-SFR bursts once every 10⁸ yr, with zero SFR in between. Each burst shifts the models to slightly redder β and lower EW(Hβ), thereby causing a slight overlap between models with different f_esc.

In this context, the worst-case scenario would be to have a single young, low-f_esc burst superimposed on an old, passively evolving population which contributes substantially to the Hβ continuum (thereby lowering EW(Hβ)) while leaving the UV slope β largely unaffected, since this could in principle mimic high f_esc in the EW(Hβ)–β diagram. In Figure 7(a), we show the evolution of such a composite population, where a 1 Myr old f_esc = 0 model has been superimposed on a 700 Myr old passively evolving population. As the mass ratio between the old population and the short burst is increased to ≳ 30, this f_esc = 0 model can wander into the region normally occupied by leaking galaxies with f_esc ≳ 0.5. If the mass ratio is increased to >500, the passively evolving population becomes completely dominant also in the UV and shifts β to very red values. However, the presence of an old, massive component gives rise to other tell-tale spectral signatures at the wavelengths probed by both NIRSpec (0.6–5 μm) and the MId-Infrared Instrument (MIRI; 5–27 μm in imaging mode) onboard the JWST.

Zoom In Zoom Out Reset image size

This is demonstrated in Figure 7(b), where we plot the SEDs of a 10 Myr old, constant-SFR model with f_esc = 0.5 (blue line) and a 1 Myr old f_esc = 0 model superimposed on a passively evolving population with age 700 Myr (black line), where the old component is 300 times more massive than the young population. Both spectra have been redshifted to z = 7. These two models occupy a very similar position in the EW(Hβ)–β diagram (black triangle and orange 10⁷ yr age marker in Figure 7(a)), but the continuum flux observed at ≳ 3 μm would be very different in the two cases (much higher in the case of the two-component galaxy).

The continuum flux at ⩽5 μm comes for free when measuring β and EW(Hβ) with JWST/NIRSpec and JWST/MIRI can trace the continuum to even longer wavelengths. MIRI lacks the sensitivity to provide useful spectroscopic limits on the continuum, but MIRI imaging should yield competitive constraints. This is demonstrated in Figure 7(b), where we include the expected detection limits for JWST/NIRSpec low-resolution (R = 100), 1–5 μm spectroscopy (red dashed line) and compare these with the corresponding limits for JWST/MIRI imaging (red triangles) at ≈5–10 μm. In both cases, we adopt an exposure time of 10 hr and require a signal-to-noise ratio (S/N) = 10 detection.

In the quest to derive f_esc from the rest-frame UV/optical SED, dust attenuation is potentially a more serious problem than metallicity or star formation history, since this part of the SED may contain insufficient information to gauge the impact of LyC extinction. Many recent studies of reionization-epoch galaxies suggest that these objects have suffered very little dust attenuation (A(V) ≲ 0.2 mag; e.g., Finkelstein et al. 2012b; Bouwens et al. 2012b; Dunlop et al. 2013; Wilkins et al. 2013; but see Schaerer & de Barros 2010 for a different view). On the other hand, it does not necessarily take much attenuation in the rest-frame optical or non-ionizing UV to have a considerable impact on the LyC—it depends on the spatial distribution of the dust compared with the gas and stars.

3.7.1. Dust Distributions

Predicting the effects of dust on galaxy SEDs is notoriously difficult, since the effective dust attenuation law depends not only on the amount and composition of dust, but also on the spatial distribution of stars, dust, and gas (e.g., Calzetti et al. 1994; Gordon et al. 1997; Witt & Gordon 2000). In Figure 8, we schematically illustrate four different dust distributions (two for each of the gas geometries of Figure 1). This is not meant to be an exhaustive list of the possible dust distributions, but serves to illustrate a few different scenarios. In Figures 8(a) and (b), the dust is distributed in a uniform screen outside the H ii region. In the cases depicted in Figures 8(a) and (c) (ionization-bounded nebulae with holes), dust has no direct impact on the leaking LyC flux, since the leakage is assumed to take place through holes that are devoid of both gas and dust. In the density-bounded nebula of Figure 8(b), on the other hand, the dust screen has a very pronounced effect on the escaping LyC. As an example of this, consider a dust screen giving a rest-frame V-band extinction of A(V) = 0.2 mag. Extinction curves are notoriously uncertain in the far-UV, but if one adopts the average LMC attenuation curve presented by Cartledge et al. (2005), this A(V) converts into a predicted LyC depletion factor of ≈3 at 910 Å.⁴ Hence, the LyC escape fraction would effectively be limited to f_esc < 1/3.

Zoom In Zoom Out Reset image size

In the case of a foreground dust screen (Figures 8(a) and (b)), the rest-frame optical attenuation can straightforwardly be estimated from JWST/NIRSpec spectra of z ≈ 6–9 galaxies by measuring of the ratios of Balmer emission lines (e.g., Hβ/Hγ) since the intrinsic ratio can be estimated from recombination theory (Hβ/Hγ ≈ 2.1–2.2 for a wide range of electron densities and temperatures, assuming Case B recombination). However, if the dust is spatially mixed with the H ii gas as in Figures 8(c) and (d), some fraction of the LyC photons will be absorbed by dust before they have a chance to ionize the gas (e.g., Inoue et al. 2001; Inoue 2001) and this will not be revealed by the Balmer line ratios. While radiation pressure and dust sublimation act to produce a dust cavity around the ionizing stars, neither of these effects is sufficiently strong to prevent LyC photons from getting absorbed directly by dust in local H ii regions (Inoue 2002). Indeed, Hirashita et al. (2003) estimate that ≈40% of the LyC photons produced in low-redshift starburst galaxies are typically absorbed by dust.

3.7.2. Dust Effects in the EW(Hβ)–β Diagram

Since simulations of high-redshift galaxies suggest that when LyC escapes it does so through essentially dust-free channels (Gnedin et al. 2008; Razoumov & Sommer-Larsen 2010), we will limit the discussion to radiation-bounded nebulae with holes (Figures 8(a) and (c)). For both of these geometries, the effects of small amounts of optical attenuation (A(V) = 0.2 mag) on positions of galaxies in the EW(Hβ)–β diagram are illustrated in Figure 9, assuming the Cartledge et al. (2005) average LMC attenuation law to capture the transport of radiation through dusty regions.

Zoom In Zoom Out Reset image size

At low f_esc, a foreground dust screen (attenuation vector represented by the solid arrow) has a pronounced effect on β but almost no effect on EW(Hβ), since the stellar continuum emerging through the screen still dominates the rest-frame UV, while the unattenuated continuum that escapes through the holes remains subdominant at wavelengths close to Hβ. At high f_esc, the reverse is true: β is largely unaffected by dust since the unattenuated stellar UV continuum escaping through the holes dominates over the dust-attenuated continuum and EW(Hβ) is lowered due to the increasing importance of stellar Hβ continuum emerging through the holes. Hence, the attenuation vector changes direction and amplitude depending on where in the EW(Hβ)–β diagram one started off.

If the dust is spatially mixed with the ionized gas (Figure 8(c)), the stellar non-ionizing continuum will still suffer A(V) = 0.2 mag of attenuation, but the Balmer line ratios will reflect a value smaller than this, since part of the attenuation happened prior to photoionization. At fixed f_esc, the nebula will also appear fainter due to LyC extinction. Here, we model this case by assuming that the star light entering the H ii region is attenuated by A(V) = 0.1 mag before it has a chance to ionize the gas (lowering the number of LyC photons that enter the gas by a factor of ≈2) and apply a further A(V) = 0.1 mag correction outside the H ii region (giving a total of A(V) = 0.2 mag, as before). Other ratios of stellar to nebular attenuation can certainly be considered, but this serves to exemplify the EW(Hβ) and β trends. The resulting LyC depletion factor is somewhat lower than suggested by Inoue (2001) and Inoue et al. (2001) for local H ii regions, but similar to that derived by Hirashita et al. (2003) for low-redshift starburst galaxies. Our procedure gives rise to the dashed attenuation arrows in Figure 9. As in the case of the outer dust screen, the β slope is reddened more at low f_esc than at high f_esc, but the overall reduction in nebular emission (due to LyC extinction) increases the relative impact of stellar Hβ continuum on the emerging SED and decreases EW(Hβ) substantially at all f_esc.

At fixed A(V), the two types of attenuation vectors discussed above have slightly different directions and amplitudes. Moreover, a measurement of the Hβ/Hγ Balmer decrement will tend to underestimate the overall attenuation whenever LyC photons are lost due to dust extinction, as in the case where ionized gas and dust are mixed. For a given Hβ/Hγ ratio, there will therefore be lingering uncertainties in how to correct the position in the EW(Hβ)–β diagram unless the amount of LyC extinction can be assessed. For example, consider an object observed to lie close to the 10⁷ yr marker on the orange line (f_esc = 0.5) in Figure 9. To correct the position of this object for dust effects, one needs to pick the right attenuation vector and move the object in the opposite direction in the diagram. If the Balmer line ratios suggest A(V) = 0.2 mag, a correction assuming an outer dust screen would place the object on the green line (implying f_esc = 0.7), whereas a correction assuming mixed gas and dust would still place it on the orange line (i.e., f_esc = 0.5). These uncertainties are not easy to overcome from observations in the rest-frame UV/optical alone, but supporting observations in the rest-frame mid/far-IR could improve the situation (as further discussed in Section 3.7.3).

The possibility that the effective attenuation law may have a more pronounced effect on emission lines than on the far-UV/optical continuum represents yet another complication. This phenomenon is observed in local starburst galaxies (e.g., Calzetti et al. 1994, 2000) and is most readily interpreted as due to temporal evolution in the opacity of star-forming regions (e.g., Charlot & Fall 2000)—the youngest star clusters (that contribute most of the LyC photons and the nebular emission) tend to be more strongly affected by dust than the slightly older star clusters, which contribute substantially to the far-UV/optical stellar continuum. In principle, this effect could mimic LyC leakage in the EW(Hβ)–β diagram by causing a substantial reduction in EW(Hβ) but only a modest reddening of β. However, we find that the Calzetti et al. (2000) attenuation scheme—in which the nebular emission experiences a dust optical depth a factor of ≈2 times than that of the stars—results in an attenuation vector no more extreme than those already included in Figure 9. In order for a f_esc = 0 galaxy to be misclassified as an f_esc ⩾ 0.5 object, one requires both a higher ratio of nebular to stellar dust opacity—as in some of the more extreme objects discussed by Wild et al. (2011)—and substantial nebular attenuation (A(V) > 0.3 mag, as revealed by, e.g., the Hβ/Hγ emission line ratio).

3.7.3. Mid/far-IR Observations

The part of the SED required for the measurement of β and EW(Hβ) falls within the wavelength range of JWST/NIRSpec for z ≈ 4–9. The use of higher-order Balmer lines instead of Hβ (Hγ, Hδ, etc.) would extend the range to z ≳ 9, but at the expense of losing the [O iii]λ5007 line, which is useful for getting a handle on the metallicity (Section 3.5). Since the proposed method is applicable for relatively blue SEDs (β ≲ −2.3) and positive EW(Hβ), the most challenging aspect of these observations is to measure the Hβ continuum. Using the JWST/NIRSpec exposure time calculator version P1.5 with the R = 100 setting, we estimate that one should be able to get a S/N ∼ 4 measurement of the Hβ continuum (giving an error σ(log₁₀EW(Hβ)) ≈ 0.1, since the error on the line flux will be negligible by comparison) by integrating over ≲ 50 Å bins in the rest-frame spectrum in 10 hr for a 100 Myr old, Z = 0.020, f_esc = 0.5, constant-SFR galaxy with stellar mass 3 × 10⁸ M_☉ at z = 7. At z = 9, the corresponding mass limit would be 1 × 10⁹ M_☉.

Galaxies this massive are well within the detection limits of ultradeep Hubble Space Telescope imaging and many ∼10⁸–10⁹ M_☉ objects have already been discovered in current z ≳ 7 samples (e.g., Schaerer & de Barros 2010; McLure et al. 2011). Strong lensing by foreground galaxy clusters may boost the fluxes of high-redshift background objects by magnification factors of up to μ ≈ 30–100 (e.g., Johansson et al. 2012; Zackrisson et al. 2012; Coe et al. 2013). By targeting lensed fields, it should therefore be possible to push these detection limits down to ≈1 × 10⁷ M_☉ (z = 7) and ≈3 × 10⁷ (z = 9), if μ ≈ 30 is assumed. In terms of the rest-frame 1500 Å luminosities often used when measuring the luminosity function at high redshifts, this corresponds to M₁₅₀₀ ≲ −16.0 at z = 7 and M₁₅₀₀ ≲ −17.5 at z = 9. These detection limits convert into SFRs of ≳ 0.1 M_☉ yr⁻¹, placing them slightly above the level where stochastic sampling of the stellar IMF starts to have serious effects on the LyC flux of galaxies (Forero-Romero & Dijkstra 2013).

On a somewhat longer timescale, similar observations will also be possible with ground-based telescopes like the Giant Magellan Telescope,⁵ the Thirty Meter Telescope,⁶ and the European Extremely Large Telescope.⁷

As discussed in Section 3.7, complementary observations of the rest-frame mid/far-IR with SPICA or ALMA could in principle provide information on the fraction of LyC photons destroyed directly by dust. Such observations could also reveal the presence of heavily embedded components within the target galaxies. Formally, the f_esc derived from JWST observations can only reflect the star formation that actually contributes to the rest-frame UV/optical SED and extremely dust-obscured star formation (undetectable in the optical, but bright at IR wavelengths, as in IR-luminous galaxies at lower redshifts; e.g., Choi et al. 2006) would shift this f_esc estimate away from that relevant for simulations of galaxy-driven reionization, where the total SFRs of galaxies need to be converted into a LyC flux emitted into the IGM. Based on the dust SEDs for star-forming dwarf galaxies presented by Takeuchi et al. (2005), we estimate that strongly lensed (μ ≈ 100), z ≈ 6–9 galaxies with SFR ∼10 M_☉ yr⁻¹ should be detectable with ALMA in ∼10 hr exposures at ≈400–1000 μm. SPICA may be able to detect such lensed z ≈ 6 galaxies at ≈80 μm in less than 1 hr, but this requires that one can push the photometry a factor ≈5 below the formal SPICA confusion limit—for instance by using auxiliary, high-resolution data from other wavelength bands (e.g., obtained with JWST or ALMA) to subtract off the flux contribution from nearby interlopers.

In its current form, the proposed method can be used to assess f_esc for the brightest galaxies at z ≳ 6, but is likely to provide interesting limits for galaxies exhibiting very high LyC escape fractions only (f_esc ≳ 0.5). However, the rest-frame UV/optical SED carries more information than contained in the EW(Hβ)–β diagnostic and stronger constraints can possibly be set by making use of the entire SED obtained from JWST/NIRSpec observations. In future papers, we intend to explore the true f_esc limits that can be obtained from spectroscopy using more detailed SED simulations. While it is likely that f_esc estimates for less extreme leakers may be obtained this way, the luminosity limits of the method are unlikely to change significantly. Hence, only objects with M₁₅₀₀ ≲ −16.0 (stellar masses ≳ 10⁷ M_☉) can be probed this way.

Galaxies in this mass range may be insufficient to explain the ionization state of the universe at z ≲ 9 (e.g., Finkelstein et al. 2012a; Ferrara & Loeb 2013; Alvarez et al. 2012; Paardekooper et al. 2013) and methods based on photometry rather than spectroscopy (in the vein of Ono et al. 2010; Bergvall et al. 2013; Pirzkal et al. 2012, 2013) will be necessary to estimate f_esc for objects at lower masses. While photometric methods are unlikely to provide strong constraints for individual targets, they can, on the other hand, be applied to larger samples of galaxies. In the future, we therefore also intend to investigate the likely f_esc limits that can be set using photometry with the JWST/Near Infrared Camera and JWST/MIRI filters.

Population III galaxies may exist at high redshifts (e.g., Stiavelli & Trenti 2010) and are predicted to have high f_esc (Johnson et al. 2009; Benson et al. 2013). Such objects may be detectable with JWST in strongly lensed fields (Zackrisson et al. 2012) and should be identifiable because of their unusual colors (Inoue 2011; Zackrisson et al. 2011a, 2011b) and spectra (e.g., Schaerer 2002; Inoue 2010, 2011; Zackrisson et al. 2011b), at least for ages up to ∼10⁷ yr (Zackrisson et al. 2011b). While the low masses predicted for such galaxies (total stellar mass perhaps no more than ∼10⁴ M_☉; Safranek-Shrader et al. 2012) makes it unlikely that these objects would contribute substantially to the reionization of the universe, the spectral diagnostics we propose could in principle be used to assess f_esc for these galaxies as well. This is demonstrated in Figure 10, where we show the evolution predicted in the EW(Hβ)–β diagram for a population III galaxy experiencing a brief burst (10⁷ yr) of zero-metallicity star formation, in the case of a radiation-bounded nebula with holes. The characteristic masses of population III star formation are generically predicted to be top-heavy, but the exact stellar IMF remains highly uncertain. Here, we have used the Raiter et al. (2010) stellar SEDs for a log-normal IMF with characteristic mass 10 M_☉, dispersion 1 M_☉, and wings extending from 1–500 M_☉, but the results are similar for other IMFs explored by Zackrisson et al. (2011b) during the early phase (∼10⁷ yr) when population III galaxies should be uniquely identifiable. At low ages, the strong impact of nebular emission on the UV slope at these low metallicities makes the model tracks curve upward in the EW(Hβ)–β diagram, but the tracks corresponding to different LyC escape fractions still remain clearly separated, thereby allowing an estimate of f_esc. At older ages (≳ 10⁷ yr), EW(Hβ) drops while β evolves very little, which causes substantial degeneracies between f_esc and the IMF (not shown, to avoid cluttering).

Zoom In Zoom Out Reset image size

One of the most important problems in observational cosmology is identifying the sources that reionized the universe. It is generally assumed that a substantial fraction of the ionizing photons were produced by galaxies, but other objects such as quasars (e.g., Dijkstra et al. 2004), mini-quasars (e.g., Madau et al. 2004), micro-quasars (e.g., Mirabel et al. 2011), and population III stars (e.g., Venkatesan et al. 2003) may also have played a role. Currently, there exist mostly indirect constraints on the timing of the reionization epoch. The spectra of high-redshift quasars suggest that reionization was completed by z ≈ 6 (Fan et al. 2006; Mortlock et al. 2011), while observations of the cosmic microwave background can be used to put limits on the extent of the reionization period (e.g., Komatsu et al. 2011; Larson et al. 2011; Pandolfi et al. 2011; Zahn et al. 2012). The detailed reionization history at z ≳ 6.5 remains highly uncertain, with statistics of Lyα emitters seemingly suggesting that the IGM was still highly neutral at z ≈ 6.5–7 (e.g., Ouchi et al. 2010; Pentericci et al. 2011; Jensen et al. 2013), while IGM temperature measurements instead indicate that reionization was already mostly complete by these redshifts (Theuns et al. 2002; Raskutti et al. 2012).

The most promising prospect for constraining the reionization history in the near future is by observing the 21 cm line emission from the partly neutral IGM. Radio observatories such as LOw Frequency ARray (LOFAR; Van Haarlem et al. 2013), the Murchison Widefield Array (MWA; Tingay et al. 2013), and Precision Array for Probing the Epoch of Reionization (PAPER; Parsons et al. 2010) are aiming for a detection at some point within the next few years. The 21 cm signal encodes a wealth of information about the physics of the reionization epoch, but the first generation of experiments will be focused on simple statistics such as the variance and power spectrum of the 21 cm brightness temperature. These statistics will not be able to constrain the evolution of the ionized fraction directly, but will have to be interpreted by comparing with simulations (e.g., McQuinn et al. 2007; Lidz et al. 2008; Iliev et al. 2012).

Most simulations to date include only galaxies as sources of ionizing photons and use simplistic recipes to assign ionizing fluxes to the sources, tuned to reproduce the available observational constraints. With measurements of the LyC escape fraction in combination with observed UV luminosity functions, much of the "wiggle room" for these recipes would be eliminated. This, in turn, would make 21 cm measurements more powerful for constraining the sources of reionization by enabling better comparisons between the 21 cm statistics predicted by different source models.