HETDEX Public Source Catalog 1—Stacking 50,000 Lyman Alpha Emitters^∗

Though the sky subtraction methodology is described in detail in Gebhardt et al. (2021), because of its direct relevance to several topics in this work, we present a brief overview of the salient points. The goal of the sky subtraction is to remove background light while preserving the light from astronomical sources. The background light comes from the atmosphere (literally, the "sky"), as well as instrumental effects such as thermal noise and scattered light in the optics. Furthermore, any large diffuse astrophysical source will be included in the sky background; for example, light from unresolved faint foreground and background galaxies that are uniformly distributed in the sky will add to the estimate of the background. There are additional complications discussed in Gebhardt et al. (2021) such as wavelength distortion in the spectrographs that can lead to sky background residuals.

Here we make use only of the local form of the HETDEX sky subtraction, as it is more stable than the full field sky subtraction in the pipeline version used to create HPSC-1, the source of the spectra for this work. The on-sky size of the local sky is 50'' × 12'', whereas the full sky is derived from the entire 21' diameter focal plane. While both sky subtractions use similar methods, the local sky subtraction employs 112 fibers for each amplifier (4 per IFU), and the full field sky subtraction uses all fibers in an exposure (∼35,000). The full sky necessitates additional adjustments for amplifier-to-amplifier variations and other complications. For emission-line source detection, we rely on the local sky since it is more robust to instrumental effects and generates smaller background residuals.

For the local sky subtraction, fibers containing obvious continua are removed, including those with counts more than 3× the biweight scale of the fibers on the amplifier. The remaining fibers are then used to compute the per wavelength bin background from their biweight locations. Typically, 70%–80% of the fibers in an amplifier are used for this calculation, though it can be as low as 25%–30% when there is a star or other bright object in or near the field of view of the amplifier. As HETDEX detections can include fibers from more than one amplifier, due to their size and/or their position within an IFU, the fibers included in their detection aperture can have different local sky subtraction corrections. One potential consequence of the local sky subtraction, which we revisit later, is that it can partly remove low-level flux actually associated with detected galaxies if that flux is spatially extended over a large portion of an IFU.

Similar to the background residual correction presented in Davis et al. (2021), the LAE spectra stacked in this work have an observed-frame residual correction applied prior to their shift to the rest frame. The purpose of this correction is to account for small systematics in the calibration, and the average extra light in the extraction apertures that come from faint, undetected background and foreground sources that are undetectable at the individual detection level. In Davis et al. (2021), apparently empty individual sky fibers are selected, sorted, and stacked on a per-observation basis and used to correct the individual fibers of the LAE detection from the matched observation. However, in this work, given the large number of statistics of the 250× larger LAE sample, we select entire apertures rather than individual fibers and average over many observations instead of applying separate corrections for each observation.

Using the observations from 2018 June, when the IFU array became significantly populated, through June 2020, the end of data acquisition for HPSC-1, we collect up to 200 spectra from random empty 35 radius apertures within each of the ∼2000 observations in this range; this radius matches the size of the standard HETDEX aperture (Gebhardt et al. 2021), and the spectra are extracted using the same method as the normal HETDEX detections. The center of each aperture is randomly selected from the footprint of the observation's IFU coverage with the following constraints:

• 1.  
  No aperture center can be within 15 of another. This means the apertures can overlap but not near their central regions where the PSF weights are highest.
• 2.  
  There must be at least 15 fibers included in the aperture over the three exposures comprising the observation. This avoids apertures where a significant portion falls off the edge of an IFU.
• 3.  
  The measured g magnitude in the PSF-weighted spectra extracted under the aperture must be fainter than 24, and have a median flux density between 3900 and 5400 Å greater than −4 × 10⁻¹⁸ erg s⁻¹ cm⁻² Å⁻¹. This satisfies the empty requirement, given that the HETDEX depth is close to 25 in g in good conditions but can be near 24 in poor seeing (Gebhardt et al. 2021; Davis et al. 2023; Mentuch Cooper et al. 2023). The negative lower limit accounts for the possibility of some over-subtraction and for excursions due to noise.
• 4.  
  There can be no detected emission lines within the extracted spectrum and no HETDEX detections within 20.

These criteria yield more than 300,000 PSF-weighted spectra (not all observations yield the target 200 empty apertures) that are stacked along their native HETDEX 2 Å wavelength bins using the biweight measure of central location (Beers et al. 1990) using the same method described in Davis et al. 2021. The result is shown in Figure 1 with the residual spectrum flux in blue and the flux uncertainties for each wavelength bin in gray. The reported uncertainties are statistical and are defined as ${\sigma }_{b}/\sqrt{N}$, where σ_b is the biweight scale of the N contributing empty aperture measured spectra for that wavelength bin. The spectrum is also made available in electronic form. Sky subtraction residual spectra created by stacking under alternate subselections of the full ∼300,000 sample are stable with respect to variations in the seeing FHWM, throughput and instrument response, date, and observed decl.

Zoom In Zoom Out Reset image size

Each LAE spectrum is corrected by subtracting this sky subtraction residual in its observed frame. For individual spectra the correction is less than 1%, far below the noise level; the systematic becomes meaningful only when stacking hundreds or thousands of galaxy spectra. Figure 1 shows a distinct rise in the far blue caused by instrumental limitations and is slightly negative redward of about 4000 Å, indicating a small over-subtraction by the reduction pipeline. Both these effects are more than 100× smaller than the typical HETDEX flux limits (∼2 × 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹, but can rise up to 10× higher blueward of 3800 Å; Gebhardt et al. 2021) and are inconsequential for individual spectra. They only become significant when stacking large numbers of spectra. Future data releases will include improvements in the data-reduction pipeline to better address these issues.

As described in the previous subsection, the spectra are first corrected for the sky subtraction residual in their observed frame. Depending on the needs of the analysis, all spectra are corrected for wavelength- and redshift-dependent fractional IGM transmission. The IGM correction applied in this work is from CIGALE (Boquien et al. 2019), which is derived from Meiksin 2006. The individual application of the IGM transmission correction can be highly uncertain since IGM attenuation is dependent on the particular sight line (Inoue et al. 2014). However, since we are averaging over many thousands of spectra, this becomes naturally statistical (Dijkstra et al. 2007; Zheng et al. 2010; Laursen et al. 2011; Behrens et al. 2018; Steidel et al. 2018; Gurung-López et al. 2019; Byrohl & Gronke 2020). Lastly, the spectra are shifted to their own rest frame so that they can be aligned for stacking.

The rest-frame shift has two components, one for the wavelength bins and one for the flux density at each bin center. Since the reported HETDEX wavelengths are in air, the observed wavelengths are converted to vacuum (Greisen et al. 2006) and then simply adjusted by (1+z) where z is determined from the Gaussian-fitted line center of the Lyα emission line (Mentuch Cooper et al. 2023). As discussed in Davis et al. 2021 and later in this work (Section 4.3), there is no correction applied for any Lyα velocity offset from an individual galaxy's systemic redshift. Depending on the science needs, this wavelength-only shift, by itself, may be sufficient, but it only represents the observed flux at rest wavelengths: when stacking, this can create a bias against higher redshift objects due to increased cosmological dimming. For this work, and to counteract this bias, we also convert the observed flux density (f_λ ) to a luminosity density (L_λ or L_ν , as needed). Please note that here we use luminosity density as an analog to flux density, not as a luminosity per unit volume. This conversion is a straightforward application of

$$\begin{eqnarray}&&{L}_{\lambda }=4\pi \,{D}_{L}^{2}\,{f}_{\lambda }(1+z),\end{eqnarray} \tag{ 1 }$$

where D_L is the luminosity distance. Where L_ν is needed, we simply multiply L_λ by λ²/c. Unless stated otherwise, z (specifically, 1 + z) is defined as the ratio of the observed frame, vacuum-shifted, Gaussian-fitted emission-line center wavelength to the vacuum rest-frame wavelength.

Once the previously described corrections and shifts are made, the spectra are stacked using the same method described in Davis et al. (2021). The rest-frame wavelength spacing is adopted from the highest redshift object to be stacked, with the blue and red wavelength endpoints coming from the highest and lowest redshift objects, respectively. All spectra are then linearly interpolated onto this wavelength grid and stacked in each wavelength bin using the weighted biweight statistic (Davis et al. 2021), as a modification of the biweight statistic described in Beers et al. (1990). The weighted biweight alters the biweight location slightly by including an additional weight for each element in the sample. The weight used in the stack is the inverse of the uncertainty on the flux measures within each wavelength bin, such that the more uncertain fluxes have a reduced contribution to the biweight location. We note that for these large data sets (≫10⁴ spectra), the biweight and weighted biweight perform very similarly to a median average, with a brief comparison of the two methods provided in Section 4. In all cases, the default tuning constants are used. For the biweight location (or weighted biweight location), the constant is 6 and for the biweight scale it is 9.

While it is common practice to normalize the spectra before stacking, often to the flux near 1500 Å for LBGs and LAEs, (Vargas et al. 2014; Marchi et al. 2017; Steidel et al. 2018, and others), the nature of the HETDEX data makes this impractical. For the vast majority of the spectra presented in this work, only the Lyα emission is detectable above the flux limits, so there is no measured continuum against which to normalize. Instead, we stack the un-normalized spectra first, and then, where appropriate for the analysis, we normalize the stack against the flux over a section of its wavelengths.

Since the HETDEX wavelength window is fixed at 3500–5500 Å and the LAEs span 1.9 < z < 3.5 (Gebhardt et al. 2021), only the wavelength region immediately around rest-frame Lyα can receive contributions from all spectra (the number of contributing spectra for each wavelength bin is included as a column in the stacked spectrum data file). The fluxes in wavelength bins increasingly blueward of Lyα are populated by LAEs at increasingly higher redshifts, while the opposite is true moving redward of Lyα.

Before proceeding further, it is necessary to acknowledge several caveats that can affect the analyses and interpretations.

Contamination of the sample of HPSC-1 LAEs by misclassified emission, in particular that of [O ii], is quite low at an estimated value of ∼2.5%–3.0% (Mentuch Cooper et al. 2023). However, since it is not zero, it may slightly dilute the signal of the underlying continuum. Another 0.5%–2.0% of contamination could come from other individually detected lines such as C iii], C iv, or Mg ii, some of which are associated with possible AGNs. And though we make efforts to exclude AGNs through the identification of broad-line (>1200 km s⁻¹) pairs (e.g., Lyα + C iv, C iv+ C iii], O vi+ Lyα, etc.), emission-line profile fitting, catalog matching, and visual inspection (Liu et al. 2022; Davis et al. 2023), some AGNs (particularly single narrow-line Type II AGNs and some broad-line objects where a second line is not identified by the pipeline) may remain in the LAE sample (Liu et al. 2022; Lujan Niemeyer et al. 2022; Davis et al. 2023; Mentuch Cooper et al. 2023). Some fraction of the detections, particularly at lower emission-line S/N, may also be false detections of noise. While work is ongoing to better quantify the specific rates of these false detections, for the LAEs of HPSC-1, as described in Section 6.6 of Mentuch Cooper et al. (2023) positively confirm that 91% of the detections in a subsample of LAEs with repeat observations and thus sets an upper limit of 9% on the total fraction of contaminants.

Unlike Davis et al. (2021), the spectra in the stack have not been selected to exclude LAEs with nearby neighbors, as that information is absent from this data release. Similarly, the spectra have not been deblended to remove the contribution of flux by neighboring objects (D. Davis et al. 2023, in preparation). While the removal of the sky subtraction residual should largely handle the average contribution of faint line-of-sight interlopers, excess flux from these sky-adjacent neighbors can still be present in the individual galaxy spectra and find its way into the stacks.

While the HETDEX catalog has no galaxy preselection, its emission-line detections are flux limited (∼4 × 10⁻¹⁷ erg s⁻¹ cm⁻²; Gebhardt et al. 2021), and the HPSC-1 catalog excludes detections with emission-line S/N < 5.5 (Mentuch Cooper et al. 2023). There is thus a bias against detecting galaxies with observationally fainter Lyα, which necessarily increases with increasing redshift. However, the use of the median-like weighted biweight mechanics (Beers et al. 1990; Davis et al. 2021) in the stack does mitigate the influence of bright outliers, which helps maintain the representative nature of the stack.

As previously stated and further discussed in Section 4.3, the spectra are aligned for stacking based on their Gaussian-fitted Lyα line. This can lead to a ∼1 Å smearing of spectral features due to the individual Lyα redshift offsets from the galaxies' systemic redshifts.

Lastly, as noted in Section 3.4, the wavelength regions farther from Lyα have fewer contributing spectra and thus increased uncertainty in the stack compared to wavelengths closer to the Lyα line. Furthermore, we have to assume that evolution in the galaxies over those redshift ranges of each stack is minor, in order to consider the stacked spectrum, as a whole, as representative of the underlying sample.

Table 1 presents a summary of several properties (see the table note for descriptions) of the full data set and full stack along with three additional subselections that divide the sample by redshift. The redshift distribution, with mean z = 2.6, for this work is shown in Figure 5. Additional binnings based on other properties will be presented in future works. With the exception of the average r magnitude, which is included from various archival photometric imaging catalogs that overlap the HETDEX observations (Davis et al. 2021, 2023), the properties are measured from the individual spectra or the stacked spectra. A comparison of the biweight location ($\tilde{{x}_{b}}$) values of g, Lyα Luminosity, and Lyα equivalent width to the values derived from the stacked spectra demonstrates that they agree very well, suggesting the stacks are, indeed, similar to an average representation of the samples. We acknowledge that both the biweights of g and W_λ (Lyα) are more limited comparisons as HETDEX cannot measure continuum magnitudes fainter than ∼25 in g; hence, W_λ (Lyα) is a lower limit (Gebhardt et al. 2021; Davis et al. 2023; Mentuch Cooper et al. 2023).

Zoom In Zoom Out Reset image size

Table 1. Summary of Ensemble Properties by Redshift Range

Redshift	N	〈g〉	Stack g	〈r〉	〈LLyα 〉	Stack LLyα	〈Wλ (Lyα)〉	Stack Wλ (Lyα)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
1.9 < z < 3.5	52,000	$\tilde{{x}_{b}}$ = 25.9, σb = 1.30	25.6	$\tilde{{x}_{b}}$ = 25.5, σb = 1.16	$\tilde{{x}_{b}}$ = 42.92, σb = 0.212	42.83	$\tilde{{x}_{b}}$ = 93.1, σb = 84.0	114.0
2.0 < z < 2.5	21,000	$\tilde{{x}_{b}}$ = 25.8, σb = 1.35	25.3	$\tilde{{x}_{b}}$ = 25.3, σb = 1.21	$\tilde{{x}_{b}}$ = 42.86, σb = 0.211	42.77	$\tilde{{x}_{b}}$ = 89.9, σb = 86.5	108.8
2.5 < z < 3.0	18,000	$\tilde{{x}_{b}}$ = 25.8, σb = 1.27	25.6	$\tilde{{x}_{b}}$ = 25.4, σb = 1.12	$\tilde{{x}_{b}}$ = 42.91, σb = 0.187	42.84	$\tilde{{x}_{b}}$ = 80.9, σb = 72.6	101.7
3.0 < z < 3.5	11,000	$\tilde{{x}_{b}}$ = 26.2, σb = 1.18	26.2	$\tilde{{x}_{b}}$ = 26.0, σb = 0.83	$\tilde{{x}_{b}}$ = 43.03, σb = 0.180	42.96	$\tilde{{x}_{b}}$ = 107, σb = 81.4	130.0

Notes. (1) Redshift range of (sub)selection. (2) Number of galaxies in the (sub)selection to the nearest 1000. (3) The biweight location ($\tilde{{x}_{b}}$) and biweight scale (σ_b ) of the Sloan Digital Sky Survey (SDSS) g magnitude of individual LAE detections as computed from the HETDEX spectra using the speclite Python package (Kirkby et al. 2023; see also, Davis et al. 2023). Though always computed, values ≳25 are fainter than the HETDEX detection limit. (4) The SDSS g magnitude computed from the observed frame stack of HETDEX spectra, again using speclite. Given the high S/N of the stack, the formal error on the fit magnitude is <0.01. (5) The biweight location and biweight scale of the r magnitudes of individual LAE detections as computed from photometric imaging with r coverage. Depth is catalog dependent (Davis et al. 2023), with 75% of HPSC-1 ≳ 26. Nondetections in the imaging are included as their respective limits. (6) The log₁₀ of the biweight location and biweight scale of the luminosity (erg per second) of the Lyα line of the individual LAE detections. (7) The log₁₀ of the luminosity (in erg per second) of the Gaussian-fitted Lyα line of the stack. Given the high S/N, uncertainties on the fits are ≪1%. (8) The biweight location and biweight scale of the rest-frame equivalent width (Å) of Lyα in the individual LAE detections. Since the continuum is often undetected, this is a lower limit. (9) The rest-frame equivalent width (Å) of Lyα in the stack. Given the high S/N, the errors on the fit are <0.05 Å.

Download table as:  ASCIITypeset image

Since we align the individual spectra based on their fitted, rest-frame Lyα emission-line centers, our figures show the stacked Lyα line centered at our adopted wavelength of 1215.67 Å. However, due to the complexities of Lyα radiative transfer and the suppression of the flux near and just blueward of rest-frame Lyα (Verhamme et al. 2006; Smith et al. 2018; Verhamme et al. 2018; Byrohl & Gronke 2020, and many others), we are often fitting to the red peak of Lyα, and there is an offset with respect to each galaxy's systemic redshift. In the stack, this manifests as a slight offset between the expected and observed positions of the other emission and absorption features. Additionally, as the Lyα offset from systemic is variable by galaxy, there can also be a smearing/broadening of these other (stacked) spectral lines. With only Lyα detected in the vast majority of the HETDEX LAE spectra, and generally at S/N ≲ 6 for HPSC-1, our ability to correct for the velocity offset of individual galaxies prior to stacking is limited (Davis et al. 2021). However, as we estimate the typical velocity offset to be at most a few hundred kilometers per second (see below), the impact on this work is small and we do not refine it further here.

To estimate the average Lyα velocity offset in our ensemble, we repeatedly fit the center lines of several emission and absorption features using a simple Markov Chain Monte Carlo (MCMC) Gaussian fitting code and compute a velocity offset from the assumed fiducial wavelength. These features are selected as they are clear, have high S/N, and are not significantly blended with any other lines (or can have the blended portion easily masked). The results are summarized in Table 2.

The velocity offsets from the features are similar, though there is some obvious scatter likely as a combination of interstellar medium (ISM) and IGM confusion, outflows, and other kinematics. The overall mean, 235 ± 18 km s⁻¹, provides a good estimate of the Lyα velocity offset for stack and thus for the typical HPSC-1 LAE. Where stated for the remainder of this work, we adopt a rounded value of 250 km s⁻¹ for the velocity offset from systemic for the stack. Though LAEs can exhibit a sizeable difference in individual Lyα velocity offsets, this average is consistent with those of LAEs and LBGs found in Shapley et al. (2003, ∼360 km s⁻¹), Erb et al. (2014, ∼240 km s⁻¹), Shibuya et al. (2014, ∼230 km s⁻¹), Steidel et al. (2018, ∼300 km s⁻¹), and Muzahid et al. (2020, ∼170 km s⁻¹). A more rigorous investigation is presented in L. Weiss et al. (2023, in preparation).

These estimates generally agree with and bracket the ∼200 km s⁻¹ reported for the much smaller, z > 3 sample in Davis et al. (2021). This equates to a less than 1 Å offset in the adopted Lyman continuum region, 880–910 Å. That work finds no significant impact on the Lyman continuum estimate measured over 30 Å and no apparent change when applying a correction such as that in Verhamme et al. (2018), Byrohl et al. (2019), and Gurung-López et al. (2020). This small velocity offset is also consistent with the expectation of enhanced Lyman continuum leakage (Izotov et al. 2018; Naidu et al. 2021).

On either side of the Lyα emission line, there are deep, negative troughs. While absorption is expected near the Lyα emission, the depths of the troughs, as shown in Figure 2, are enhanced by our reductions. Though scattering of photons near Lyα resonance by the H i in and around the LAEs is expected, the majority of the depths of these troughs are a result of the HETDEX sky subtraction and reduction pipeline (Section 3.1). The sky subtraction assumes that the sky off-source is the same as the sky on-source. Any differences in this assumption could then create a feature in the resulting stacked spectra. As one possible scenario, the version of the pipeline in this data release could over-subtract faint Lyα emission from halos around the LAEs that extend ∼10 or more arcseconds. In this case, we could create a self-subtraction of the Lyα emission and cause the troughs to go negative. Another possible scenario is that a diffuse UV background exists as part of the general sky background we measure. In this case, the on-source sky could be different than the off-source sky due to absorption around Lyα. A more complete discussion of the artificial troughs and real Lyα scattering is presented in L. Weiss et al. (2023, in preparation). For the measurements of this work, we simply avoid the Lyα troughs by masking them.

The calculations of the Lyα luminosity for both individual galaxy spectra and the stack are similar. For an individual galaxy, the identified Lyα emission line is fit with a simple Gaussian (Davis et al. 2023) whose area is the integrated flux, noting that the continuum level is a free parameter and allowed to be negative. That flux is then converted to a luminosity using Equation (1), but without the tailing (1 + z) term. For the L_λ stacked spectra, since we are already in the rest frame and in luminosity density, we simply fit a Gaussian where the continuum is set by the L_λ redward of Lyα and whose area is the integrated line luminosity with the Lyα trough regions (Section 4.4) masked out. This is the observed luminosity of the escaping Lyα. There are no additional corrections applied for extinction or attenuation by the rest-frame dust and ISM.

As shown in Table 1 and Figure 6, the individual galaxy Lyα log luminosities (i.e., log₁₀(L_Lyα /[erg s⁻¹])) in this sample range from 42.26–44.45 with a biweight location of 42.92 ± 0.212 and 1.9% of the Lyα emission lines greater than 43.50. At the brighter end, it is probable that some of the galaxies host unidentified AGNs and some luminosities may be inflated due to measurement uncertainty. The full sample-stacked spectra have a Lyα luminosity of 42.83, about 20% lower than the sample biweight, though well within the uncertainty. The fit to the Lyα line of the stack is a more precise average of the sample as a whole due to the increased S/N, since the continuum level is a free parameter for both the stack and the individual detections. The stacked continuum is well detected where continuum is rarely detected in individual HETDEX LAEs and the greatest spread in luminosities in the sample is associated with the lowest Lyα S/N detections, which represent the largest fraction of the sample (Figure 6).

Zoom In Zoom Out Reset image size

Table 1 naively shows a weak trend of increasing Lyα luminosity with increasing redshift bin, which could indicate some small evolution with redshift (Ciardullo et al. 2011, for example). However, the decrease is consistent with the loss due to the HETDEX flux limits and cosmological dimming.

The Lyα rest-frame equivalent widths, W_λ (Lyα), are not explicitly reported in HPSC-1 but are shown in Figure 7, along with the sample biweight location and the corresponding stack value. Each W_λ (Lyα) is computed from the integrated Lyα flux and a combined estimation of the continuum from the spectrum and associated photometry (Davis et al. 2023; Mentuch Cooper et al. 2023) divided by (1 + z). For each of the rest-frame stacked spectra (Table 1), W_λ (Lyα) is simply the ratio of the integrated Lyα luminosity to the fitted continuum around the Lyα line, again with the Lyα troughs masked. As with the Lyα luminosities, the stacked spectra measurement of W_λ (Lyα) is consistent with the biweight measure from the sample, but the stack may provide a more robust average as the continuum is detected; HETDEX rarely detects the continua of z ≳ 1.9 galaxies.

Zoom In Zoom Out Reset image size

The HPSC-1 equivalent width distribution is similar to that found with MUSE in Kerutt et al. 2022. HPSC-1 reports that 15% of its LAEs have W_λ (Lyα) >240 Å and the full sample has a biweight location value of 93.1 ± 84.0 Å; Kerutt et al. 2022 reports that 16% (in their full sample) of the galaxies have W_λ (Lyα) >240 Å and the characteristic W_λ (Lyα) is 95.5 Å. Restricting the comparison of the two samples to their LAEs within the overlapping redshift coverage, 2.9 < z < 3.5, the distributions deviate a bit more, but remain similar, see Figure 8. The redshift-restricted MUSE sample of 591 LAEs (out of 1920 LAEs of the full sample) has 16% with W_λ (Lyα) above 240 Å versus 11% for the redshift-restricted HPSC-1 sample of 13,500 LAEs. The biweight locations of the restricted MUSE W_λ (Lyα) and HPSC-1 samples are 85.4 and 102.7 Å, respectively.

Zoom In Zoom Out Reset image size

The HPSC-1 equivalent widths given in Table 1 and Figure 9 might also suggest some weak evolution with redshift similar to findings in Ciardullo et al. 2011. However, this is even less clear with the HPSC-1 data than the possible trend with Lyα luminosity as the equivalent width measures are less secure due to the lack of detected continua and are also biased toward the detection of brighter Lyα at higher redshifts due to cosmological dimming.

Zoom In Zoom Out Reset image size

The observed UV continuum slope, β, modeled as f(λ) ∝ λ^β , is taken from the full stack (Figure 2) and measured from 1250–1850 Å by fitting a simple least squares optimized power law after masking the features near 1260, 1302, 1335, 1394, 1549, 1640, and 1665 Å. Given the wavelength range, this is similar to β₁₈ in Calzetti 2001. The best-fit β, −2.36 ± 0.09, implies that our average LAE has a young, low-metallicity stellar population, a large ionizing photon production efficiency, and very little dust extinction (Calzetti 2001; Bouwens et al. 2010b; Calabrò et al. 2021; Chisholm et al. 2022; Saldana-Lopez et al. 2022). This is in line with that expected for the z ∼ 3 LAE population (see also Section 4.9). This also suggests an increased likelihood of Lyman continuum escape and may indicate a recent (5–15 Myr) burst of star formation (Calzetti 2001; Calabrò et al. 2021; Chisholm et al. 2022; Saldana-Lopez et al. 2022).

The absolute UV magnitude (M_UV or M₁₅₀₀) of the full stack is computed from the mean and median luminosity density between 1400 and 1600 Å in the rest frame, using the standard 3631 Jy AB magnitude zero-point scaled to L_ν,0 of 4.3454 × 10²⁰ erg s⁻¹ Hz⁻¹ at 10 pc. The M_UV from the mean values corresponds to a medium bright galaxy with −19.67 ± 0.13, while the median values are essentially identical, −19.68 ± 0.13. The errors are statistical and come from the propagation of the standard deviation of the flux density between 1400 and 1600 Å in the rest frame.

P Cygni line profiles (Beals 1934) are clearly visible in Figure 2, for example O vi (1032,1038 Å), N v (1241 Å), and, to a lesser degree, C iv (1549 Å). These profiles are complex combinations of nebular and stellar absorption and emission, and contain a wealth of information on the stellar population, star formation history, and the metal enrichment of the galaxy, but a proper decomposition and fitting is beyond the intended scope of this work. However, qualitatively speaking, the P Cygni profile is a telltale indicator of strong stellar winds, massive stars, and recent star formation (Steidel et al. 1996; Leitherer et al. 2001; Chisholm et al. 2019). While not exhibiting a P Cygni, the extremely broad He ii (1640 Å) emission, 940 ± 130 km s⁻¹, is also strongly suggestive of a very young population (Chisholm et al. 2019). This is consistent with our picture of LAEs, and the clarity of the profile in the stack further highlights the S/N gains and additional physics that is not directly accessible in the individual HPSC-1 spectra.

To model and gain some understanding of the underlying stellar population, we perform spectral energy distribution (SED) fitting on the stacked spectrum (Figure 2). Since we are limited by the wavelength coverage in this spectrum to the rest-frame UV, we are only probing more recent star formation. The SED fitting is performed with the Python package for Fitting the stellar Continuum of UV Spectra or FiCUS ²⁵ (Saldana-Lopez et al. 2022). Four separate runs use the Starburst99 (SB99; Leitherer et al. 2010) and the Binary Population and Spectral Synthesis (BPASS, v2.2.1; Eldridge et al. 2017) single-burst stellar population models, along with the dust attenuation models from (Reddy et al. 2016, hereafter R16) and the SMC extinction model by Prevot et al. (1984). All runs assume the same Kroupa (2001) initial mass function (IMF) with a 100 M_⊙ upper limit and fit for 10 stellar ages (1, 2, 3, 4, 5, 8, 10, 15, 20, and 40 Myr) and four metallicities (5%, 20%, 40%, and 100% of Z_⊙). We use the IGM corrected version of the full stack spectrum (Figure 2) with an applied velocity offset of 250 km s⁻¹ (Section 4.3) and fit over 915–1915 Å in the rest frame, assuming a mean redshift of 2.604 for the 50,000 contributing galaxies. Key results of the four runs are summarized in Table 3 with the fit with the best χ² (row 2, SB99 + SMC) shown in Figure 10. We are not modeling Lyman continuum escape with FiCUS, as we are only fitting the continuum redward of the Lyman limit and are not yet confident using the restricted wavelength range for the z > 3.0 HETDEX LAEs where Lyman continuum could be observed.

Zoom In Zoom Out Reset image size

While there are some small differences in the results from the four runs, with BPASS favoring slightly older, less enriched stellar populations and the R16 dust model favoring slightly increased reddening, a consistent picture emerges of an ensemble-averaged galaxy with a young (10–15 Myr), metal-poor (0.2–0.3 Z_⊙) stellar population with minimal (E(B − V) 0.03–0.10) reddening. The observed and intrinsic UV continuum slopes (ranging from β = −2.10 to −2.05 and −2.65 to −2.54, respectively) for the four runs are fit over 1250–1850 Å as described in Section 4.7 and bracket the −2.36 UV continuum slope measured directly from the HPSC-1 stack. Saldana-Lopez et al. 2022, with support from Izotov et al. (2016), Salim et al. (2018), and Shivaei et al. (2020), suggest that SB99 + SMC may best represent the conditions for these LAEs and, indeed, though perhaps coincidentally, that combination does produce the lowest χ² fit. This model (SB99 + SMC) does favor the youngest, UV light-weighted average stellar age with a significant fraction, just over 75%, of the light coming from stars with ages less than 5 Myr.

For this measurement, we use L_ν instead of f_ν as explained in Section 3.4. We subselect the 11,000 HPSC-1 LAEs with z > 3, where the rest-frame Lyman continuum region, defined as (880–910 Å) for consistency with the literature (Shapley et al. 2006; Marchi et al. 2017; Steidel et al. 2018, and many others), fallts within the HETDEX spectral range. The L₉₀₀ value used is the weighted biweight location of the 68 wavelength bin luminosities in the residual subtracted (Section 3.2) and IGM transmission corrected (Boquien et al. 2019) L_ν spectrum between 880 and 910 Å. The reported error is statistical; a standard error analog as the biweight scale of the same data divided by $\sqrt{n}$, where n = 68.

Since 1500 Å is not available in the HETDEX spectral window for z > 3 galaxies, the L₁₅₀₀ value is estimated using two methods. The first estimate is based on a scaling of the weighted biweight flux density between 1268 and 1296 Å using the slope fit between 1250 and 1525 Å from our 2.6 < z < 3.5 stack of 24,000 LAEs. That wider redshift range is selected as the narrowest window that includes sufficient wavelength coverage to the red of Lyα to capture the L₁₅₀₀ region. The L₁₅₀₀ normalized fitted slope is 5.11 (± 0.561) × 10⁻⁴ L_ν Å⁻¹, and results in a stack average IGM corrected L₉₀₀/L₁₅₀₀(out) = 0.169 ± 0.0361 (or L₉₀₀/L₁₅₀₀(obs) = 0.069 ± 0.0137, without the IGM correction) as an upper limit, given the aforementioned caveats. The labels, out and obs, adopt the convention in Steidel et al. (2018) and Pahl et al. (2021, 2023). The second L₁₅₀₀ estimate simply uses the weighted biweight luminosity density of the 113 wavelength bins between 1475 and 1525 Å taken directly from the 2.6 < z < 3.5 stack using the same method as L₉₀₀ described earlier. This result differs by ≪1%, completely consistent with the first estimate. Repeating the same computation, but substituting the median and standard deviation for the weighted biweight location and biweight scale, yields an IGM transmission corrected L₉₀₀/L₁₅₀₀(out) = 0.159 ± 0.0455, or L₉₀₀/L₁₅₀₀(obs) = 0.066 ± 0.0172 when not correcting for the IGM.

To place this measurement in some context, we compare against the stack of the subset (26 out of 124 galaxies) of the LBG-selected, z ∼ 3 galaxies in the Keck Lyman Continuum Spectroscopic Survey (KLCS; Steidel et al. 2018) that are also classical LAEs with (W_λ (Lyα) > 20 Å (Pahl et al. 2021). The overplotted stacks are presented in Figure 11. Both stacks are normalized to their own flux near 1500 Å and interpolated onto the same wavelength grid and smoothed with a 1 pixel (0.44 Å) Gaussian kernel (σ). The HETDEX (HPSC-1) stack has also been shifted to correct for its approximate 250 km s⁻¹ Lyα velocity offset (Section 4.3). Though the individual spectra contributing to the KLCS stack are of higher S/N with longer exposures, they have similar resolving powers, with R ∼ 800 for HETDEX and for KLCS at λ < 5000 Å and R ∼ 1400 for KLCS at λ > 5000 Å (Steidel et al. 2018; Gebhardt et al. 2021). Our L₉₀₀/L₁₅₀₀ estimates are defined in a similar way as the <f₉₀₀/f₁₅₀₀>_out measurements of Steidel et al. 2018; Pahl et al. 2021 and others.

Zoom In Zoom Out Reset image size

As a brief aside, the Lyα troughs (Section 4.4) shown in the HPSC-1 stack are much deeper than in the KLCS stack. The KLCS LAE subsample is comprised of UV-bright, LBG-selected galaxies with weaker Lyα, as compared to the majority of the 50,000 HPSC-1 LAEs. This difference in selection may probe galaxies with somewhat different halo and internal properties which could contribute to some differences in the troughs. However, as discussed in L. Weiss et al. (2023, in preparation) and briefly in Section 4.4, while the existence of the Lyα troughs is physically motivated, their manifestation within the HPSC-1 stack is substantially enhanced by the HETDEX data-reduction pipeline.

While our relative Lyman continuum values, L₉₀₀/L₁₅₀₀(out), are ∼1.5–2.0× higher than what is found in Pahl et al. (2021) and Pahl et al. (2023) for the galaxies in their highest (W_λ (Lyα) > 20 Å) equivalent width bin, based on a comparison to the properties of the HPSC-1 sample, we might expect a larger escape of Lyman continuum photons. Additionally, the z ≥ 3 for ∼900 Å rest-frame IGM transmission in the model used in this work (Boquien et al. 2019) is about 10% higher than that in Steidel et al. (2018) and Pahl et al. (2021), making the IGM correction here slightly smaller. As can be seen in Figure 11, the equivalent width of the HPSC-1 Lyα line is roughly 3× larger than that of the KLCS LAE subsample, with HPSC-1 measuring 114 Å (or 130 Å for the 3.0 < z < 3.5 subsample, see Table 1) and the KLCS LAE subsample measuring 36 Å using the same MCMC fitting code used in Section 4.6. Measuring the UV continuum slope as in Section 4.7, we find the KLCS subsample stack with a shallower, but still very blue, β = −1.88 versus the −2.36 of the HPSC-1 stack. The FiCUS analyses for the KLCS subsample stack, with the same configurations as in Section 4.9, suggest a similarly young stellar population component with similar metallicity but a somewhat larger E(B – V) as compared to the HPSC-1 sample. Using the SB99+SMC configuration, though all 4 runs show the same relative differences, we see the UV light-weighted average stellar age = 9.36 ± 1.39 Myr, Z(Z_⊙) = 0.21 ± 0.02, and E(B – V) = 0.052 ± 0.001.

Both the KLCS and the HPSC-1 LAEs show consistent ages for their recent star formation as well as comparable metallicities, suggesting similar ionizing photon production. However, the stronger Lyα emission, lower E(B – V), and steeper UVβ slope in the HPSC-1 sample promote the expectation of increased leakage of those photons from the HPSC-1 LAEs (Behrens et al. 2014; Verhamme et al. 2015; Smith et al. 2018; Calabrò et al. 2021; Naidu et al. 2021; Chisholm et al. 2022; Maji et al. 2022; Saldana-Lopez et al. 2022).

We emphasize that, while this relative L₉₀₀/L₁₅₀₀ measurement can be used in context with other works, it is only a step in obtaining an estimate of the average, intrinsic escape fraction of ionizing photons from these galaxies. We also caution again (see Section 4.1) that the HPSC-1 sample is not as strictly controlled for contamination from on-sky neighbors as in Pahl et al. (2021) and that may influence the results. The subtraction of the average empty aperture (Section 3.2) helps compensate for contributions of light from faint interlopers, but there is no exclusion of LAEs with detected, sky-adjacent neighbors as in Davis et al. (2021) nor is there any deblending of light from such neighbors as is performed in D. Davis et al. (2023, in preparation). As such, this result is an upper limit. A more complete and extended study of the Lyman continuum escape from z > 3 LAEs will be presented in D. Davis et al. (2023, in preparation).