Detection of Lyman Continuum from 3.0 < z < 3.5 Galaxies in the HETDEX Survey

• 1.  
  Emission Wavelength: Because the spectral region of interest is in the Lyman continuum (880–910 Å) and the HETDEX observable wavelength range is 3470–5540 Å, the minimum observed emission-line wavelength is selected as 4860 Å (or z ≈ 3 for Lyα; see panel 2 in Figures 1 and 4). This ensures the blue end of our targeted Lyman continuum range is captured by the VIRUS spectrographs, allowing for some small variation in each detector's wavelength range and avoiding the edgemost pixels. The wavelength limits for the Lyman continuum are chosen to select the photons responsible for reionization (those blueward of the Lyman limit) in the region surrounding the host galaxies without excessive contamination from outside sources. The (red) upper bound is fixed by the Lyman limit itself while the (blue) lower limit is a bit softer and is established by the mean free path of Lyman continuum photons in 3 ≤ z ≤ 4. The optical depth of the IGM at z ∼ 3 results in a mean free path of Δz ≈ 0.18, which in turn sets the blue limit to ∼870 Å at z = 3 and ∼880 Å by z = 4 (Haardt & Madau 1998; Rudie et al. 2012). Because the 880–910 Å range is commonly used in other publications, it is adopted here as well for consistency.
• 2.  
  Emission-line FWHM and g-band Apparent Magnitude: To exclude AGN and brighter contaminants, the maximum allowed fit FWHM to the emission line is set at 600 km s⁻¹ with a bright g-band magnitude limit of 23. While AGN in the redshift range of interest are certainly LAEs, in this work we are interested in the far more common, simple star-forming galaxies (though AGN and broad-line/bright LAEs will be included in future work). This selection criterion does not absolutely guarantee the complete exclusion of AGN (particularly faint, narrow-line Type II AGN), but the 600 km s⁻¹ limit is well below the more typical ∼1000–2000 km s⁻¹ used in classifying Type I (broad-line) AGN (Antonucci 1993; Urry & Padovani 1995; Steidel et al. 2002; Baron et al. 2016, and others). Future work will explore the possible contamination by Type II AGN, but here the impact is assumed to be negligible. The photometric g coverage for this sample is incomplete and the g magnitude used here is computed by passing the HETDEX collected spectrum through the SDSS g filter (Doi et al. 2010) using the Python speclite package (version 0.8; Kirkby 2020). The estimated apparent g magnitude is limited by the HETDEX flux sensitivity to 24.5–25.0 (depending on the observing conditions, exposure time, and the VIRUS IFU) (KG21), and most of the bandpass-related computation is based on additional photometric imaging (generally r for this work).
• 3.  
  S/N and χ² Model Fit: Although there is good recovery (at the time of this writing, the detailed HETDEX survey recovery rate at the lower S/N range is still being evaluated) down to an S/N of approximately 5 for the Lyα emission line, and the absolute number of LAEs grows rapidly with decreasing S/N (Figure 1 panel 1), the relative number of false detections also increases with decreasing S/N and, in an effort to keep this preliminary manual examination sample manageable, an emission-line S/N minimum threshold of 6.0 is enforced. Closely related to the S/N cut, a χ² limit of 1.2 on the single Gaussian emission-line fit is imposed to remove potentially questionable detections. While there certainly are real sources with a poor χ² fit (particularly broad, noisy, or strongly asymmetric lines) this criterion performs generally well. This selection, combined with the redshift range and emission-line FWHM cuts, resulting in a moderately bright Lyα line luminosity distribution centered near 1.2 × 10⁴³ erg s⁻¹ (panel 4 in Figure 1).
• 4.  
  P(Lyα) Classification: As an obvious condition, candidates that are not strongly indicated as LAEs are excluded from the initial selection with a minimum required value of 0.8 on the ELiXer composite P(Lyα) classification value. While expressed as a real number (0–1), P(Lyα) is not a proper probability but does represent a form of confidence in the classification of the emission line as Lyα. The full details of the computation are presented in DD21. In brief, however, it is a weighted combination of an analysis incorporating several EW estimates, redshifts, and line fluxes (Farrow et al. 2021), based on the Bayesian analysis in Leung (2017). Also factoring in are emission-line positions, flux ratios, and other estimated physical and spectral properties from the HETDEX spectra and available photometry. The Lyα EW is estimated as

  $$\begin{eqnarray}&&\mathrm{EW}=\displaystyle \frac{F}{C(1+{z}_{\mathrm{Ly}\alpha })},\end{eqnarray} \tag{ 3 }$$

  where F is the integrated line flux as fit by HETDEX (KG21) and C is the continuum flux density as estimated, in this work, from either the r aperture magnitude or the SDSS g magnitude from the HETDEX spectrum, as described earlier. While the Lyα line appears within g for the redshifts of this work, g-band photometry is often not available or is not as deep as r and so the r magnitude is used, assuming a flat continuum between Lyα and r, as the primary continuum estimate in the EW calculation. This does not alter the efficacy of the EW analysis on which the emission-line discrimination is based and, in fact, may produce a cleaner separation (Adams et al. 2011; Leung 2017). The Lyα EW from the HETDEX spectra is more similar to a proper EW in that it estimates the continuum nearer the emission line (as opposed to the r magnitude). The sample distribution of this EW estimate is shown in panel 3 of Figure 1 with the caution that the continuum estimates include wavelengths blueward of Lyα, subject to increased IGM attenuation, as well as the Lyα line itself. These two issues at least partly cancel but are sources of error in the EW estimates. Future work will refine the EW estimates and the current results are shown here only in the context of the broad description of the data set. With that caveat, the smallest rest-frame emission-line EW in the final data set is ∼40 Å for the r-continuum estimate and ∼26 Å for the HETDEX spectrum estimate, both greater than the commonly used 20 Å EW minimum as an LAE discriminator (Cowie & Hu 1998; Gronwall et al. 2007; Adams et al. 2011). In approximately 15% of the final sample, the imaging counterpart cannot be detected by the Python photometric package (sep, Barbary 2016). However, in all cases, the single emission line is obvious with the assumption that the counterpart is undetected because it is fainter than the imaging limit (here, 26.5 is adopted for HSC-r; Figure 1, panel 5). The classification of the emission line as Lyα based on the EW analysis, when the continuum is fainter than 25th magnitude, is supported (Leung 2017; DD21; KG21).
• 5.  
  Spatially Resolved Candidates: LAEs at z ∼ 3 are expected to be a few kiloparsecs to maybe 10 kpc in radius (Ribeiro et al. 2016) and should be generally unresolved in the HSC-r imaging (seeing FWHM 06–10). Any HSC-r detected galaxies with an sep elliptical aperture effective radius or ELiXer circular aperture radius (DD21) ≥20 (≳15 kpc for z ∼ 3) are rejected as they could be misclassified as an AGN or foreground galaxy or include a slightly offset line-of-sight interloper, extending the projected profile. Here the effective radius is defined as

  $$\begin{eqnarray}&&\displaystyle \frac{1}{2}\sqrt{{ab}},\end{eqnarray} \tag{ 4 }$$

  where a is the semimajor axis and b is the semiminor axis of the aperture ellipse. It is also possible that the candidate galaxy is actually two or more interacting LAEs (this is also true even when the candidate is unresolved), but without the higher spatial resolution of HST imaging, it is difficult to impossible to distinguish these cases. In regions of the HETDEX survey with multiband HST coverage this becomes possible using SED comparisons as in (Pahl et al. 2021), but given the large HETDEX fiber size (K. Gebhardt et al. 2021, in preparation), it would challenging to separate the individual LAE spectra. These could be common and are certainly interesting cases for reionization as interactions can lead to enhanced star formation (Keel et al. 1985; Lawrence et al. 1989; Jogee et al. 2009, and others).
• 6.  
  Neighbor Flux Contamination: Regions around known, nearby large galaxies (such as M101 in a region defined by R.A. (J2000) ∈ [2102, 2113] and decl. (J2000) ∈ [541, 548]) are completely excluded from consideration to avoid the risk of contamination from flux contributed by the outer regions of foreground galaxies. Additionally, no candidate is permitted to have any detected object (using the sep Python package) extending to within 2'' (resolved or unresolved) of the HETDEX detection position. This reduces the chances of emission contamination and source confusion, with the astrometric (centroid) accuracy ∼ 05 and the FWHM seeing generally better than 2'' (KG21).For slightly larger separations, the PSF-weighted flux of all imaging found sources within 9 arcsec² of the candidate LAE are summed and the candidate LAE is removed from the sample if the potential flux contamination exceeds a threshold. As a note pertaining to the execution sequence, this check is performed after the initial query conditions described above and prior to the manual inspections described in the next subsection. The flux for each neighboring source (and the LAE candidate) is computed with simple aperture photometry (again, in HSC-r) and convolved with the HETDEX PSF for the observation of the candidate LAE. For simplicity, each source is assumed to be unresolved (noting that larger, extended sources within this area would automatically cause the LAE candidate to be rejected) with a flat spectrum over the bandpass. The (moffat) modeled PSF-weighted fraction (KG21) of flux in the overlap with the LAE candidate and every other source is summed over and the candidate is rejected if that sum (over all on-sky neighbors within the 9 arcsec²) is larger than 25% of the candidate LAE's r flux. The choice of the 25% threshold is somewhat rough and is intended to be similar to the average uncertainty in the aperture magnitudes that provide the flux estimates. For the final sample of 214 galaxies, the biweight midvariance location of the fractional PSF-weighted flux overlap is 0.01 ± 0.015, far below the 25% limit, and that potential flux contribution is ignored. More sophisticated modeling with a flux correction may be used for future samples to allow some relaxation of this criterion.

These initial selection criteria and automated conditions reduce the set of detections from over 1.5 million to approximately 1000 for visual inspection.

The final step in producing the data set for this work is the manual inspection of the remaining detections that pass the previous selection conditions. Despite best efforts to impose rigorous, repeatable conditions, this is not free from observer subjectivity. Future data selections with improved automation (in development at the time of this writing) will seek to minimize this subjectivity.

The visual inspection examines the following conditions and excises detections that meet any of the following:

• 1.  
  Poor Imaging: If the photometric imaging within the few arcseconds around a candidate LAE is unclear, corrupt, incomplete, or shallow (magnitude limit < 25), it cannot be reliably examined for the potential foreground interlopers. Any such detections with poor imaging are removed.
• 2.  
  Data Issues: The 2D HETDEX spectra are examined for cosmic-ray strikes, hot pixels, bad columns, or other data corruption in the Lyman continuum spectral region in the nearest four fibers (in order of increasing distance from the detection center position) and in the stacked sum of all contributing fibers. Any issues in these cutouts result in the candidate's exclusion. Similar data issues in the Lyα region do not warrant removal unless they cast doubt on the actual detection or classification as Lyα.
• 3.  
  Bad Neighbors: As an extension of the Neighbor Flux Contamination condition in the previous subsection, the photometric imaging cutouts for each remaining LAE candidate are examined for the presence of neighboring sources that could be contributing unwanted flux, including any suggestion of very faint, low-surface-brightness interlopers that could be missed by the automated source extraction. Larger area (30'' × 30'') zoom-outs of the images are also examined for any potential contaminants (generally, stars or small, but spatially resolved, nearby galaxies) that could be sources of scattered-light/flux contamination. Again, any candidates with suspect neighbors are excluded.
• 4.  
  Unidentified Emission Lines: The presence of any unidentified emission lines—that is, any visually apparent or ELiXer Gaussian fitted (DD21) emission lines in the spectra that are inconsistent with the classification as an LAE at the assumed redshift or suggestive of a second, blended spectrum from a lower-z line-of-sight interloper—in the full 1D spectrum suggests some contamination (blending) of spectra. While the ELiXer application performs automated line searching and fitting, this manual check is an additional safety against any convincing emission lines that ELiXer fails to identify. Candidate LAEs with such spectra are rejected.

Following this examination, 214 detections remain to comprise the final data set for this investigation (again, see Figure 1 for summary statistics). Figure 2 presents a subselection of 15 candidate LAEs as they appear in the HSC-r imaging with Table 1 listing the basic positional and brightness data for those galaxies (with the full listing made available in a machine-readable format). The LAE candidates all appear compact and point-source like, with no apparent foreground interlopers.

Zoom In Zoom Out Reset image size

"Sky fibers" are effectively empty fibers (i.e., those that do not fall on any apparent astrophysical object), and they play an important role in this work. After HETDEX processing, any residual flux in these fibers is ostensibly due to random noise, including instrumental noise, and imperfect sky measurement and subtraction. While the majority of fibers in the typical HETDEX exposure fall on empty sky (KG21), we are interested in those freest from any measurable light from the PSF wings of nearby objects or otherwise scattered along the light path as the most representative of the true sky. Functionally, sky fibers, within each three-exposure set, are identified by the following procedure that acts on the post-sky-subtracted, rectified spectra:

• 1.  
  To remove fibers with a possible measured continuum, any fibers with spectra exhibiting an average flux density (over 500 Å wide bins) outside of ±5 × 10⁻¹⁸ (erg s⁻¹ cm⁻² Å⁻¹) or ±4 × 10⁻¹⁸ (erg s⁻¹ cm⁻² Å⁻¹) over (almost) the full width (3600–5400 Å) are eliminated. The 3600–5400 Å range is selected to avoid the edges of the detector and the prominent sky lines near 3545 and 5462 Å. The choice of bounding flux densities is an empirical definition and based on the HETDEX flux limits. Because most (50%–70%) fibers fall on empty sky, the values are selected to be small enough to exclude HETDEX-detectable continuum and potential oversubtraction of sky (KG21) and still include roughly half of all fibers for a typical HETDEX exposure set. Here, by "typical" we mean those exposure sets without large errors (instrument failures) or significant foreground objects (bright stars or relatively nearby galaxies, like M101) that cover a large fraction of the fibers.
• 2.  
  Fibers with potential emission lines are also removed. For expediency, rather than attempt to fit ∼1000 Gaussian profiles per single-fiber spectrum (with ∼100,000 spectra per exposure set), here a possible emission line is defined in a sliding window of three wavelength bins, totaling 6 Å wide, with the integrated flux over the central wavelength bin exceeding 5 × 10⁻¹⁷ (erg s⁻¹ cm⁻²) and the two immediately adjacent wavelength bins exceeding 4 × 10⁻¹⁷ (erg s⁻¹ cm⁻²) each. These emission-line flux thresholds are selected based on a typical HETDEX flux limit under good conditions (noting that the actual flux limit is a function of seeing FWHM, throughput, wavelength, CCD response, fiber, etc.) (KG21).
• 3.  
  The remaining fibers are sorted by their integrated flux over 3600–5400 Å and one-third of the fibers in the sum of the two tails of this approximately normal distribution are eliminated (i.e., keeping the central two-thirds or roughly 1σ). This additional down-selection is simply an extra precaution against the influence of outliers on this residual. Figure 3 shows a representative example of the distribution of sky fiber flux (before the tails are removed). Here we display the histogram, peaking very close to zero, of the summed fluxes in the 23,000 identified sky fibers in exposure set #18 taken on 2020 June 25.
• 4.  
  The remaining fibers (roughly 30% of all fibers or about ∼25,000 fibers per exposure set) are considered effectively empty.

If the background subtraction was perfect, the distribution would center on zero (it is very close; Figure 3) and no wavelength bins or ranges of bins should vary from zero with anything but random noise. To be clear, we use the term "background" to encompass both the true sky background and instrumental noise. Figure 4 reveals, however, there is some small average residual, particularly on the blue end of the spectrum, that must be corrected. This figure (the stacked background residual spectrum from exposure set #18 taken on 2020 June 25) is typical of HETDEX observations in the fraction of sky fibers and in the residual spectrum. Even though there are larger residuals in the blue, the peak near 3550 Å and 3 × 10⁻¹⁹ (erg s⁻¹ cm⁻² Å⁻¹) is still far below the stacked, residual-subtracted galaxy flux in the Lyman continuum region presented later in this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Stacking of the spectra in the defined sky fibers (again, within each data sample included exposure set) is performed on a per-wavelength basis, with the wavelength bins aligned on their "native" observed-frame (2 Å) boundaries, using a biweight scheme, not a sum or average. The biweight measure of the central location (Beers et al. 1990) is similar to a median and is selected for its stability and robustness against outlier influence, reducing or eliminating any bias toward the relatively few, brighter sources. The biweight implementation used in this work is slightly modified with an additional weight applied as the inverse uncertainty in the flux measurements. Hereafter, we refer to this modification of the standard biweight location (ζ_biloc) simply as the "weighted biweight:"

$$\begin{eqnarray}&&{\zeta }_{\mathrm{biloc}}=M+\displaystyle \frac{{\sum }_{{v}_{i}\gt 0}({x}_{i}-M){v}_{i}^{2}}{{\sum }_{{v}_{i}\gt 0}{v}_{i}^{2}},\end{eqnarray} \tag{ 5 }$$

where M is an initial guess (typically the sample median) and the weights, v_i , are defined as

$$\begin{eqnarray}&&{v}_{i}=1-{u}_{i}^{2}\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{u}_{i}=\displaystyle \frac{{x}_{i}-M}{c\times \mathrm{median}\left\{| {x}_{i}-M| \right\}}\end{eqnarray} \tag{ 7 }$$

where x_i are the data and c a tuning constant (commonly 6; see Beers et al. 1990).

The biweight location "shifts" the median by the average difference of points to the median (x_i − M). This average excludes outliers and assigns higher weights to points close to the median than points that are farther away. As this shift is a weighted average, we can include the additional inverse uncertainty weights w_i by modifying the weights ${v}_{i}^{2}$:

$$\begin{eqnarray}&&{v}_{i}^{2}\to {v}_{i}^{2}+{\widetilde{{w}_{i}}}^{2},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&\widetilde{{w}_{i}}=\displaystyle \frac{{w}_{i}}{2\times \mathrm{median}\left\{{w}_{i}\right\}}.\end{eqnarray} \tag{ 9 }$$

This "normalization" of the additional weights, w_i , is necessary to ensure that v_i and $\widetilde{{w}_{i}}$ have comparable values and therefore similarly large contributions to the final weights. The weighted biweight location becomes

$$\begin{eqnarray}&&{\widetilde{\zeta }}_{\mathrm{biloc}}=M+\displaystyle \frac{{\sum }_{{v}_{i}\gt 0}({x}_{i}-M)({v}_{i}^{2}+{\widetilde{w}}_{i}^{2})}{{\sum }_{{v}_{i}\gt 0}({v}_{i}^{2}+{\widetilde{w}}_{i}^{2})}.\end{eqnarray} \tag{ 10 }$$

Although great pains are taken to remove the background signal from the galaxy spectra in the primary sky subtraction steps, two additional steps are taken to further refine this correction. First, to more consistently subtract the background spectrum from the galaxy spectra, each of the LAE candidate galaxy spectra are re-extracted from the original HETDEX data using the full-field sky subtraction to replace the IFU local sky subtraction used in the original iHDR2 catalog. Here, "full-field" refers to all fibers in all IFUs for an exposure set and "local" refers only to those fibers in the single IFU where the galaxy detection is made. This approach greatly expands the number of fibers used to compute the sky and reduces the potential for oversubtraction that can result from diffuse (or stray) light in a large fraction of the fibers in a single IFU (KG21; G. Zeimann et al. 2021, in preparation).

Second, for every galaxy spectrum, the independently selected sky fibers (as previously described) for all IFUs and exposures that contributed to the calculated full-field sky background for that exposure set are collected and averaged (again, using a weighted biweight). The specific PSF model for that exposure set is then applied to the resulting residual background spectrum and subtracted from the galaxy spectrum in the observed frame. Logically, this is equivalent to subtracting this background residual from the individual galaxies fibers before they are PSF weighted and combined. An example of the background residual from one exposure set (#18 from 2020 June 25) is shown in Figure 4.

In an additional check, the region around the particularly bright sky line at 3535–3555 Å is masked from the galaxy spectra. The resulting galaxy-weighted biweight stacks and Lyman continuum region averages (not shown) exhibit no discernible differences, so no extra correction for this sky line (beyond those described above) is applied.

Generally, aside from the primary Lyα emission line, the S/N of astrophysical photons over other wavelength ranges for our z ≳ 3 LAEs is too low for detection in individual sources. We therefore employ stacking to boost the signal. Because the observed Lyα line is found at different wavelengths (4860–5500 Å), the cleaned and corrected galaxy spectra are first redshifted to their own rest frames, as determined from the Lyα line, so they can be aligned for the weighted biweight stacking. The individual redshifting produces wavelength bins of slightly different widths, from 0.44 to 0.50 Å. The largest redshift results in the narrowest bins and is adopted as the grid onto which all other rest-frame spectra are linearly interpolated. The spectra are then aligned on the Lyα bin, for now ignoring systemic Lyα velocity offsets (Verhamme et al. 2018) (discussed later), and the weighted biweight is computed for each wavelength bin. The uncertainty is reported as a standard error on the biweight location, with the square root of the biweight midvariance (or biweight scale) as σ. The relatively small uncertainty, roughly 0.05 Å–0.12 Å in the rest frame, in the Gaussian fit line center (KG21) is also ignored as it is much smaller than the systemic velocity error and, unlike the velocity error, averages to zero over the sample. As with the sky fibers, averaging is performed with the modified (weighted) biweight.

Figure 5 (top panel) presents the result of stacking the full sample of candidate LAEs. Because of their different redshifts but uniform observed-frame wavelength range, the wavelength bins in the extreme outer blue and red ends of the weighted biweight stack, well beyond the Lyman continuum and Lyα regions, include fewer data points (lower panel) and are truncated from the next figures.

Zoom In Zoom Out Reset image size

Figures 6 and 7 show the main galaxy spectra stacks for this work, with the wavelength bin data listed in Table 2. To reduce the jitter and boost S/N, the final galaxy-weighted biweight stacks are summed over bins 11 elements wide (∼5 Å, given the uniform rest-frame bin width of ∼0.45 Å). The black line represents the binned stack of the entire sample. Also shown are binned stacks of the brightest one-third (in blue) and faintest two-thirds (red) by M_UV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all three weighted biweight stacks, some continuum on either side of Lyα is clearly present (at least redward of the Lyman limit), and absorption for some of the Lyman series is apparent. Other common absorption features are also marked, with several prominently represented, which will be useful for future stellar population fitting. Of note is what appears to be a P Cygni profile of [O vi] (λ1032, 1038 Å) as a signature of the young stellar populations (Leitherer et al. 1999) in these LAEs.

As is expected simply by the selection, the brightest one-third shows increased continuum over the faintest two-thirds until near and blueward of the Lyman limit, where the bright subsample falls off rapidly. This could be consistent with an expected trend to a higher EUV escape with decreasing UV luminosity (and decreasing halo mass) (Nakajima et al. 2018; Finkelstein et al. 2019), though the faint-end limit of this sample is far too bright and the sample size far too small for any conclusions. At a minimum, the absence of a statistically significant detection of EUV emission in the M_UV brightest one-third versus the clear detection in the fainter subsamples suggests an M_UV limit where EUV escape declines rapidly. Even with the decreased HETDEX sensitivity to the blue end of the spectra range, the Lyman continuum region in the full sample weighted biweight stack of galaxy spectra from this work is clearly above and inconsistent with zero flux (Figure 7), supporting the claim of flux detection in this limited sample.

Three averaging methods are applied to the Lyman continuum region with somewhat different detection significances but with consistent results. In the first method, the weighted biweight measure of central location, hereafter "weighted biweight," of the flux density of the Lyman continuum region for each individual (rest-frame) spectrum is computed and the weighted biweight of those individual EUV flux densities is calculated. For this method, the result, 0.10 ± 0.036 μJy, has the lowest significance, in part due to the larger error from the fewer wavelength bins (∼66 per rest spectrum) but is the most straightforward as the Lyman continuum regions from the individual spectra are taken as whole inputs. It is this result that is reported as the primary finding as it is the most conservative of the three methods. In the second method, all Lyman continuum region wavelength bins from all sample spectra are placed in a single vector (∼14,000 elements) and the weighted biweight is computed. The results are consistent, 0.11 ± 0.024 μJy, with a larger significance. In the third method, the weighted biweight of the Lyman continuum region of the resulting spectrum from the full stack (described above) is computed. Again, the results are consistent, 0.12 ± 0.021 μJy with the highest significance of the three methods, but due to the weighted biweight of the stacking, each wavelength bin from individual spectra does not necessarily receive the same weighting, and so a given spectrum may contribute somewhat differently to each ∼0.5 Å bin in the final stack. However, this is only a small caveat in the interpretation and does not affect the conclusions, noting that all three methods provide very similar statistics.

This same procedure is repeated for a few subselections of the data, and the results are presented in Table 3. As noted in the previous section, the trend to a higher absolute EUV emission with decreasing M_UV persists even with this small sample size and low significance. This result might support models of reionization that favor an early start with most of the ionizing radiation budget supplied by numerous, fainter star-forming galaxies (see Finkelstein et al. (2019) for an extensive discussion). Additional data will allow for substantial refinement in both S/N and luminosity bin resolution.

Table 3. Summary of Subsample Properties with EUV (900 Å) Emission

Description	Count	$\tilde{z}$	${\tilde{M}}_{\mathrm{UV}}$	Method 1	Method 2	Method 3
			±σ	(μJy) ± SE	(μJy) ± SE	(μJy) ± SE
Full Sample	214	3.16	−20.35 ± 0.711	0.10 ± 0.036	0.11 ± 0.024	0.12 ± 0.021
Brightest one-third	72	3.15	−20.99 ± 0.273	0.01 ± 0.054	0.01 ± 0.043	0.02 ± 0.037
Faintest two-thirds	142	3.17	−20.00 ± 0.532	0.15 ± 0.045	0.15 ± 0.029	0.18 ± 0.026
Faintest one-third	72	3.18	−19.40 a ± 0.327	0.13 ± 0.066	0.12 ± 0.041	0.14 ± 0.035

Notes. Summary of basic properties of the sample subselection (note: the ̃ overscore indicates a median). The redshift distribution of each subsample are effectively the same with almost identical median and scatter over [3.00, 3.49]. All reported statistics are the weighted biweight with the error on M_UV as the subsample biweight scale. The error on the three EUV escape columns are the standard errors. Method 1 computes the flux density in the Lyman continuum region of each individual candidate LAE spectrum and then performs a weighted biweight over those flux densities. Method 2 concatenates all Lyman continuum wavelength bins across all LAE candidates and performs a weighted biweight on the resulting single vector. Method 3 aligns the rest-frame wavelength bins of all candidate LAE spectra and stacks with a weighted biweight and then computes the flux density in the Lyman continuum region of the resulting stacked spectrum. See Section 4.4.

^a The faintest one-third ${\tilde{M}}_{\mathrm{UV}}$ may be fainter as 30 of the 72 galaxies are at or below the r imaging flux limit.

Download table as:  ASCIITypeset image

The complex resonant scattering of Lyα photons (Field 1959; Moses 1980, and others) means that the emission line typically separates into a blue peak and a red peak, with the blue peak often highly suppressed, resulting in a (sometimes large) line dispersion from the many scatterings and a systemic velocity offset to the galaxy when fitting only the red peak. This can create a problem when aligning the spectra for stacking and certainly for any precise characterization of the Lyα line. These rest-frame velocity offsets can be in the hundreds of km s⁻¹ (Shapley et al. 2003; Berry et al. 2012; Verhamme et al. 2018) and are generally correlated with the Lyα line FWHM. The low spectral resolution of VIRUS (R ∼ 800) amounts to an emission-line minimum FWHM resolution of ∼375 km s⁻¹ (KG21; G. Zeimann et al. 2021, in preparation) for our Lyα lines. This is sufficient for line resolution, except for the narrowest emission lines (39 of 214 Lyα lines are below 375 km s⁻¹ FWHM). Velocity offsets to the fitted line center are measurable with greater precision (∼30 km s⁻¹, or Δz ∼ 0.0005 KG21).

With the moderate range of luminosities and line widths of the LAEs in this work and the more modest objective of simply detecting Lyman continuum emission, combined with the relatively broad spectral region over which the Lyman continuum was averaged, these complications are largely ignored. With that caveat, an approximate correction for the systemic velocity offset of the red peak of Lyα (${V}_{\mathrm{peak}}^{\mathrm{red}}$) from its host galaxy is explored with Equation (2) from Verhamme et al. (2018):

$$\begin{eqnarray}&&{V}_{\mathrm{peak}}^{\mathrm{red}}=0.9(\pm 0.14)\times {\mathrm{FWHM}}_{\mathrm{Ly}\alpha }-34(\pm 60)\ \mathrm{km}\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 11 }$$

Applying this relation to the individual spectra prior to stacking, and using the 375 km s⁻¹ minimum FWHM resolution for those lines below that limit, produces a ≲0.1% overall change in individual redshifts, around 200 km s⁻¹ averaged over the entire sample (well above the HETDEX line center precision), and results in no significant difference in the Lyman continuum flux density using any of the averaging methods. Methods 1, 2, and 3 exhibit changes of +2%, +5%, and +10%, respectively. These corrections are all far less than the uncertainty in the weighted biweight of the Lyman continuum emission and have no meaningful impact on the detection at this level. With the current small sample size and moderate scatter in the Lyα line FWHM, the spectral stack is constructed from a somewhat inhomogeneous sample, and a singular correction applied poststacking would not be appropriate. However, with expanded future data sets allowing for stacks of larger numbers of more uniform luminosities, Lyα line widths, etc,, the velocity correction, in particular, could become more important in resolving ensemble spectral features and will be revisited.

Figure 8 presents an overview of the Lyα velocity correction. The top panel shows essentially the full-width view of the full sample LAE weighted biweight stack with the corrected spectrum (red) over-plotted with the uncorrected spectrum (blue). The lower two panels highlight the Lyman continuum region (left) and the Lyα line (right), respectively. The shift in the line center is clearly visible in the lower-right panel but is still relatively small and makes no significant difference over the Lyman continuum region (lower-left panel).

Zoom In Zoom Out Reset image size

Systematics in error analysis are difficult to quantify and can have a particularly large impact with smaller data sets. In an effort to report the most accurate detection significance, the following procedure (comparing subset standard errors on the mean to the standard deviations of the means of the subsets) is used to estimate a correction for the combined systematics (statistical, measurement, and instrument).

• 1.  
  Randomly subdivide the sample into (m) nonoverlapping, nonrepeating proper subsamples of size (n) (where n is ≤ an integer of the total sample size divided by m).
• 2.  
  Create a single vector of fluxes over the target wavelength range (i.e., 880–910 Å for the Lyman continuum region) across all (n) spectra for each of the (m) subsets and take the standard error of those fluxes (${\sigma }_{n}/\sqrt{n}$). Here, n is actually the subset size × the number of wavelength bins (which varies slightly from 65 depending on the subset maximum redshift), but for compactness of notation is represented by the number of galaxies in the subset.
• 3.  
  Compute the mean of each of the (m) subsample flux vectors and take the standard deviation of those means (σ_m ).
• 4.  
  Repeat steps (1)–(3) 10,000 times with different random subset draws.
• 5.  
  Compute the mean of the σ_m and ${\sigma }_{n}/\sqrt{n}$ values over those 10,000 draws and take the ratio of those means $\left({\overline{\sigma }}_{m}/({\overline{\sigma }}_{n}/\sqrt{n})\right)$.
• 6.  
  Repeat (1)–(5) for each combination of (m, n) (down to a minimum m of 2).
• 7.  
  Fit a line through the ratios (indexed by n) and extrapolate to the full data set to find the systematic correction factor (to divide into the full data set standard error).

For this data set, the ratio (${\overline{\sigma }}_{m}$/(${\overline{\sigma }}_{n}/\sqrt{n}$)) fit reaches 1:1 near (n = 82 (for m = 2)), indicating that the sample is already sufficiently large to not need an additional systematics correction applied to the standard error.

This work assumes that there is no contamination of the measured flux resulting from misclassification of the emission line, PSF scattered light from sky-adjacent astrophysical objects, or extra contributions from line-of-sight, low-surface-brightness interlopers. While multiple steps were taken to exclude these potential contaminants from the data set (see Section 3), no additional corrections, beyond the background residual subtraction (Section 4.1 ), are applied to the individual or stacked spectra. However, the galaxy stack is examined for indications of emission lines that would suggest the presence of misclassified [O ii] (λ3727) or [O iii] (λ5007), for example, within the data set. The galaxy spectra are aligned and stacked assuming the rest-frame wavelength of each potential misclassified emission line, without a correction for any velocity offset (as would be appropriate for these nonresonant line contaminants), and checked for emission lines supporting the assumption that either of these contaminants is present. No such emission is identified. Additional contaminants, [Mg ii] (λ2795+λ2802) and [C iii] (λ1909), are also specifically considered. The [C iii] line does not appear as a redshifted single line within the observed spectral range considered by this work (4860–5540 Å) and would be expected to be accompanied by [C iv] (λ1549), which is not found. The [Mg ii] doublet (at the edge of resolution as separate lines by VIRUS; KG21) likewise does not appear by itself for any part of the wavelength range of interest and should be accompanied by [C ii] (λ2326) over the entire range and by [C iii] after [Mg ii] passes 5088 Å. The [Mg ii] line is often broad and frequently (but not always) associated with AGN, and the maximum emission-line FWHM cut of 10 Å (see Section 3) and absence of continuum in individual galaxy spectra made this an unlikely contaminant. [Mg ii] can also be narrow and appear similar to Lyα; however, for the narrow cases, no evidence of either [C ii] or [C iii] is found in the sample, nor do any of the line shapes suggest a doublet peak, favoring the bluer (2795 Å) peak as is expected for [Mg ii] (Chisholm et al. 2020).

Because the continuum-selected KLCS sample median is ∼2× brighter (in r and M_UV) than that of the galaxies in our sample, the 52 brightest galaxies from this work are chosen to stack as their median magnitude (r = 24.5, M_UV = −21) matches that of the KLCS sample. We note that the M_UV distributions of the two samples are quite different, roughly symmetric about −21 for the KLCS and strongly biased to fainter magnitudes in this work (see Figure 1, panels 5 and 6). Additionally, because the HETDEX sample is specifically selected from LAEs (generally young, blue, and compact) there may be a bias for galaxies with enhanced Lyman continuum escape. For consistency with the KLCS, a flux normalization is applied to the resulting weighted biweight stacked spectrum. Due to the lower S/N of the HETDEX data, the normalization is applied to the stack and not individually as in the KLCS case. This flux normalization of the bright galaxy stack from this work approximates the f₁₅₀₀ normalization of the KLCS sample (where the KLCS spectra are individually normalized to their flux averaged over 1475–1525 Å prior to stacking). Because the rest-frame spectral range of the samples from this work does not extend to 1500 Å, a rough translation is employed using the ratio of the normalized flux densities in the KLCS sample over the fairly flat and feature free regions near 1500 Å to that near 1300 Å (which is covered in the stack from this work) as a scaling factor (specifically, between 1468 Å and 1496 Å to that of 1268–1296 Å). The bright galaxy-weighted biweight stack is normalized by its own f₁₃₀₀ × the above ratio (i.e., 0.58 ± 0.042 μJy × 1.09). While the subsample size is below the threshold where systematics can have a small impact, as the purpose of this comparison is simply to check the overall curve shape and verify key features, no systematics correction was applied.

This comparison to the KLCS sample is shown in Figure 9. For readability and to aid in the comparison, the bright subsample (shown in black) and the KLCS sample (shown in blue, without error) are each smoothed with a Gaussian kernel (σ = 2.2 Å or 5× the 0.45 Å width of the rest-frame galaxy stack wavelength bins). Although the bright subsample of this work is of a much lower spectral resolution and S/N than the KLCS sample, the two stacks align well and reproduce much the same broad shape and absorption features, lending additional validation to the methodologies of this work. Though of limited value given the large noise in this sample of 52, the observed 〈f₉₀₀/f₁₅₀₀〉 = 0.07 ± 0.182, including an IGM correction (Haardt & Madau 1995; Inoue & Iwata 2008; Steidel et al. 2018). While consistent with 〈f₉₀₀/f₁₅₀₀〉_out = 0.057 ± 0.006 in Steidel et al. (2018), it is also consistent with zero and a nondetection of Lyman continuum.

Zoom In Zoom Out Reset image size

In an update to Steidel et al. (2018), Pahl et al. (2021) use deep HST imaging to identify several previously undetected, very faint line-of-sight interlopers contaminating the KLCS sample, leading to the removal of 2/15 of the 3σ Lyman continuum detections and another 2 from the nondetections (4/124 total). This reduces the full sample ionizing emissivity in the original KLCS findings by ≲10% (from 6.0 to 5.5 × 10²⁴ erg s⁻¹ Hz⁻¹ Mpc⁻³) and the median Lyman continuum flux density in the 3σ detections by ≲5% (from 0.044 to 0.042 μJy). Pahl et al. (2021) also perform an alternate ionizing emissivity estimation, retroactively applying it to the Steidel et al. (2018) sample, and find a decrease of almost 25% (from 7.2 to 5.5 × 10²⁴ erg s⁻¹ Hz⁻¹ Mpc⁻³).

While we do not find evidence of significant foreground contamination in our sample, it is a possibility and an issue that we will continue to explore with an expanded HETDEX LAE sample. As a simple check for this work, and using the 3%–13% (4/124 and 2/15) foreground contamination rate in the previous paragraph as a rough guide, we examine the impact of similar contamination on the Lyman continuum detection. We impose a 10% contamination rate uniformly applied across our entire sample with a worst-case scenario of exact line-of-sight colocation of foreground contaminants with LAE candidates. The contaminants are simulated as simple, flat (in f_ν ) spectra with 27 in g, assuming brighter foreground galaxies with strong emission lines would be caught by the inspections in Section 3. We run 1000 iterations, randomly selecting 10% of the sample in each run, subtract the simulated contaminant spectrum from those LAE spectra, and perform the Lyman continuum stacking analyses. We find that up to 20%–25% of the flux measured in the Lyman continuum region, across all three Methods (Section 4 and Table 3), could come from undetected contamination. As a logical check, we repeat this same test, but remove the randomly impacted LAE candidates under the assumption that the contaminants are now detected. As expected, there is no significant change versus the original weighted biweight measure of the Lyman continuum emission in any of the three Methods and, due to the slightly smaller sample size, the uncertainty in those measures increases by ∼5%. The presence of undetected contamination will have an impact on future work exploring the escaping fraction of EUV photons. Some improvements in contamination detection are anticipated, particularly where deeper, multiband imaging is available, and more sophisticated simulations will be employed to account for undetected contamination, but this does not change the conclusion of this work of the detection of EUV emission in the ensemble of LAE galaxies.

As noted in the introduction, much of the practical motivation for this work is to provide an alternative survey of z ∼ 3 Lyman continuum leakers that is not subject to the same biases and limitations of other surveys that rely on target preselection (typically using the Lyman break) and deep spectroscopic observations of those targets individually. That is not to say that the methods presented here are free from biases, just that the biases are perhaps less significant taken in the context of galaxy ensemble, and some biases, notably in identifying Lyman continuum leakers (described below) are desirable for this work.

With the notable exception of the relatively small-area MUSE IFU surveys (Drake et al. 2017; Feltre et al. 2020 and others), most other z ∼ 3 surveys rely on target preselection generally using the filter dropout to identify LBGs in the desired redshift range. This is a robust selection method but does immediately introduce a selection bias that favors galaxies with a strong Lyman break, i.e., galaxies that are UV bright and Lyman continuum faint, which may tend toward older, more massive galaxies (Gawiser et al. 2007; Finkelstein 2016; Jose et al. 2013; Santos et al. 2020, 2021). UV faint galaxies, particularly those with a stronger Lyman continuum can be missed as they would be faint in the imaging with a weaker Lyman break. Thus, LBG galaxy surveys may favor galaxies with low EUV emission. This work's emission-line technique, on the other hand, has no preselection and is only flux limited and thus can frequently include fainter galaxies (as noted versus the KLCS above). It does, however, introduce the obvious bias that only galaxies with sufficiently strong Lyα emission (KG21) are included and is blind to non-LAE (or weaker LAE) galaxies. As strong Lyα emission appears common in Lyman continuum leakers, though strong Lyα emission is not sufficient, by itself, to predict Lyman continuum leakage (Kornei et al. 2010; Bian & Fan 2020), there is a potential to favor finding galaxies that are more likely to have stronger EUV emission. This does not appear to be a strong bias as evidenced by the consistent-with-zero detection of EUV emission in this work's M_UV-brightest subsample.

From the ground, except for extreme cases, z ∼ 3 galaxies are generally unresolved point sources and a slight misalignment of spectrograph singular apertures (a slit or one fiber) with respect to the source centroid can result in an under measurement of flux. Such misalignments can be the result of mechanical positioning errors of the aperture or centering on a centroid that is not colocated with the emission to be studied, though the latter should be much less of an issue given the unresolved nature of the sources and sufficiently large apertures. Since HETDEX is a blind survey, no fiber prepositioning is involved and centroids are found empirically (KG21) from the wide multifiber coverage. With the exception of sources that are detected near the edge of an IFU, which are excluded from this work, the loss of light due to fiber placement is minimized.

Many previous preselection works stack spectra from a single instrument from data collected over one or a few consecutive (or nearly consecutive) nights and thus have generally similar observing conditions and instrument states. In contrast, though individual HETDEX detections are built from three sequential 6-minute exposures, HETDEX stacks are made from observations scattered over many days, seasons, and years with varying sky conditions and instrument states. This makes background subtraction considerably more heterogeneous for HETDEX data and, though great care is taken in the background subtraction (KG21 and this work, Section 4.2), it subjects the spectra in this work to larger uncertainties. Fortunately, with the projected very large number of spectra (∼100× those of this work), the ensemble uncertainties will be minimized.

The selection criteria (Section 3), which are intended to ensure a contamination-free sample, may also introduce some bias in this work. We are primarily cautious of two criteria. First, the high EW (and other factors) used in the P(Lyα) metric (Section 3.1, bullet 4) requirement could introduce a bias against LAEs with weaker Lyα emission and stronger continuum, and thus potentially reduced Lyman continuum leakage, in an effort to avoid contamination from a misidentification of [O ii]. Second, the restriction to somewhat isolated galaxies (Section 3.1 bullet 6 and Section 3.2) could also introduce an environmental bias, though that is unclear as there is no attempt to establish the physical distances between our candidates and their on-sky neighbors. In future data sets, as improvements are made in the classification of HETDEX galaxies (DD21), the bias from the P(Lyα) selection will be reduced and added spectral deblending capabilities will help eliminate unwanted angular/spatial environmental selection effects.