# Stellar Populations of Ly$\alpha$ Emitting Galaxies in the HETDEX Survey I:
An Analysis of LAEs in the GOODS-N Field

Adam P. McCarron Department of Astronomy, The University of Texas at Austin,
2515 Speedway, Austin, TX 78712, USA<EMAIL_ADDRESS>Steven L.
Finkelstein Department of Astronomy, The University of Texas at Austin, 2515
Speedway, Austin, TX 78712, USA  Oscar A. Chavez Ortiz Department of
Astronomy, The University of Texas at Austin, 2515 Speedway, Austin, TX 78712,
USA  Dustin Davis Department of Astronomy, The University of Texas at Austin,
2515 Speedway, Austin, TX 78712, USA  Erin Mentuch Cooper Department of
Astronomy, The University of Texas at Austin, 2515 Speedway, Austin, TX 78712,
USA  McDonald Observatory, University of Texas at Austin, 2515 Speedway,
Austin, TX 78712, USA  Intae Jung Department of Physics, The Catholic
University of America, Washington, DC 20064, USA  Astrophysics Science
Division, Goddard Space Flight Center, Greenbelt, MD 20771, USA  Center for
Research and Exploration in Space Science and Technology, NASA/GSFC,
Greenbelt, MD 20771  Delaney R. White Department of Astronomy, The University
of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA  Gene C. K. Leung
Department of Astronomy, The University of Texas at Austin, 2515 Speedway,
Austin, TX 78712, USA  Karl Gebhardt Department of Astronomy, The University
of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA  Viviana Acquaviva
Physics Department, NYC College of Technology, 300 Jay Street, Brooklyn, NY
11201, USA  Center for Computational Astrophysics, Flatiron Institute, New
York, NY 10010, USA  William P. Bowman Department of Astronomy & Astrophysics,
The Pennsylvania State University, University Park, PA 16802, USA  Institute
for Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA  Robin Ciardullo Department of Astronomy & Astrophysics,
The Pennsylvania State University, University Park, PA 16802, USA  Institute
for Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA  Eric Gawiser Department of Physics and Astronomy,
Rutgers, The State University, Piscataway, NJ 08854, USA  Caryl Gronwall
Department of Astronomy & Astrophysics, The Pennsylvania State University,
University Park, PA 16802, USA  Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA  Gary J. Hill
McDonald Observatory, University of Texas at Austin, 2515 Speedway, Austin, TX
78712, USA  Department of Astronomy, The University of Texas at Austin, 2515
Speedway, Austin, TX 78712, USA  Wolfram Kollatschny  Institut für
Astrophysik, Universität Göttingen, Friedrich-Hund Platz 1, D-37077 Göttingen,
Germany  Martin Landriau Lawrence Berkeley National Laboratory, 1 Cyclotron
Road, Berkeley, CA 94720, USA  Chenxu Liu Department of Astronomy, The
University of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA  Daniel N.
Mock Department of Physics, Florida State University, Tallahassee, Florida
32306  Ariel G. Sánchez Max-Planck-Institut für extraterrestrische Physik,
Postfach 1312, Giessenbachstr., 85748 Garching, Germany

###### Abstract

We present the results of a stellar-population analysis of 72 Ly$\alpha$
emitting galaxies (LAEs) in GOODS-N at $1.9<z<3.5$ spectroscopically
identified by the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX). We
provide a method for connecting emission-line detections from the blind
spectroscopic survey to imaging counterparts, a crucial tool needed as HETDEX
builds a massive database of $\sim$ 1 million Ly$\alpha$ detections. Using
photometric data spanning as many as 11 filters covering $0.4<\lambda\
(\mathrm{\mu m})<4.5$ from the Hubble and Spitzer Space Telescopes, we study
the objects’ global properties and explore which properties impact the
strength of Ly$\alpha$ emission. We measure a median stellar mass of
$0.8^{+2.9}_{-0.5}\times 10^{9}\ \mathrm{M_{\odot}}$ and conclude that the
physical properties of HETDEX spectroscopically-selected LAEs are comparable
to LAEs selected by previous deep narrow band studies. We find that stellar
mass and star formation rate correlate strongly with the Ly$\alpha$ equivalent
width. We then use a known sample of $z>7$ LAEs to perform a proto-study of
predicting Ly$\alpha$ emission from galaxies in the Epoch of Reionization,
finding agreement at the $1\sigma$ level between prediction and observation
for the majority of strong emitters.

††software: get_spectrum.py
(https://github.com/HETDEX/hetdex_api), emcee (Foreman-Mackey et al. 2013),
Bagpipes (Carnall et al. 2018), Numpy (Harris et al. 2020)

##  1 Introduction

Lyman-alpha emitting galaxies (hereafter LAEs) have fascinated astronomers for
decades, from when Partridge & Peebles (1967) first predicted that primitive
galaxies in formation could emit a detectable Ly$\alpha$ line, through their
discovery by Cowie & Hu (1998) and Rhoads et al. (2000). These objects exhibit
strong emission of the Ly$\alpha$ photon corresponding to the $n\\!=\\!2$ to
$n\\!=\\!1$ resonant transition in hydrogen atoms. These photons are expected
to face high optical depths from neutral hydrogen to escape the galaxies in
which they are generated, and dust grains along their paths can absorb them.
To date, despite enormous effort (see Ouchi et al. 2020 for review), the
community has not formed a strong consensus on exactly how Ly$\alpha$
radiation escapes its host galaxy, and no reliable model exists to predict the
Ly$\alpha$ luminosity or equivalent width, $W_{\lambda}({\mathrm{Ly\alpha}})$,
of a galaxy given its global physical properties, such as stellar mass,
metallicity, age, star formation rate, and dust extinction.

Part of the problem arises from discrepant conclusions drawn from studying
LAEs identified using different selection techniques. Locally ($z\ll 1$), the
ultraviolet (UV) flux measured in wide or narrow-band filters often defines
LAE samples, biasing studies to brighter, higher mass systems than those found
spectroscopically (Hayes et al., 2014). Observations in the nearby universe
paint LAEs as low mass galaxies with young stellar ages as determined from
spectral energy distribution (SED) fitting, and many studies concur on trends
showing an increase in Ly$\alpha$ luminosity with decreasing dust and metals
(Hayes, 2015). Nonetheless, many galaxies show stronger Ly$\alpha$ emission
than models would predict based on dust extinction (e.g., Martin et al. 2015,
Atek et al. 2014, Scarlata et al. 2009, Finkelstein et al. 2009), and a
satisfactory explanation of this Ly$\alpha$ enhancement does not currently
exist.

With narrow-band selected LAEs at higher redshift, discrepant results still
persist. Finkelstein et al. (2009) found LAEs at $z\sim 4.5$ represent a
diverse population in terms of stellar age, mass, and dust extinction.
Finkelstein et al. (2015) modeled the SEDs of IRAC-detected LAEs at $z\sim 5$
and found a third have old stellar populations, contrasting with the young
populations found in the local universe, and Guaita et al. (2011) observed
similar heterogeneous populations in a narrow-band selected sample at $z\simeq
2.1$. Moreover, Gawiser et al. (2007) found NB-selected LAEs at $z=3.1$ to
generally be low mass, dust-free objects, but their model allowed for both
young and more evolved stellar populations, and Acquaviva et al. (2012) found
LAEs at $z=3.1$ to be older than those at $z=2.1$. Kornei et al. (2010)
compiled a UV continuum selected sample of $z\sim 3$ galaxies, finding those
with strong Ly$\alpha$ emission had older stellar populations with lower star
formation rates and less dust. Recently, Santos et al. (2020) used SED fitting
of nearly 4000 LAEs in the COSMOS field at $2<z<6$ to find that LAEs were
younger and/or more dust-poor than other UV-selected objects based on their UV
slopes.

Studies of LAE samples compiled using detection of the Ly$\alpha$ emission
line itself in the high-redshift universe confound consensus as well. Hagen et
al. (2016) used the Hobby Eberly Telescope Dark Energy Experiment (HETDEX)
pilot survey (Adams et al. 2011, Blanc et al. 2011) to compare properties of
LAEs at $z\sim 2$ with optical emission line-selected galaxies (oELGs) and
found no significant differences between the populations. Remarkably, even the
UV-slope did not differ in the two samples, implying either that diffuse dust
in the interstellar medium (ISM) did not modulate Ly$\alpha$ emission or that
oELGs strongly emit Ly$\alpha$. Recently, spectroscopic surveys have also
yielded confusing results about LAEs at $z>2$. Using data from the VANDELS
survey, Marchi et al. (2019) suggested LAEs have low mass and low dust
extinction, but found no correlation with star formation rate. From the VIMOS
Ultra-Deep Survey, Hathi et al. (2016) concurred with LAEs having lower mass
and lower dust extinction, but they found that the objects have lower SFRs
than non-LAEs. Approaching the problem from the other direction, Oyarzún et
al. (2017) found from studying the spectra of stellar mass selected galaxies
at $3<z<4.6$ that a negative correlation existed between Ly$\alpha$ equivalent
width and both stellar mass and star formation rate. A review of the field’s
current knowledge of high-redshift Ly$\alpha$ emission can be found in Ouchi
et al. (2020).

A deeper understanding of what makes LAEs unique from other star forming
galaxies (SFGs) tantalizes astronomers because of the profound implications
for leveraging LAEs as sensitive probes of reionization at $z\gtrsim 6$.
Whether the Universe re-ionized rapidly at late times (e.g. Robertson et al.
2015) or gradually beginning very early in its history (e.g. Finkelstein et
al. 2019) can determine if massive, rare galaxies or low-mass, ubiquitous
objects emitted the needed ionizing photons. Answering such a fundamental
cosmological question hinges on our ability to detect neutral hydrogen in the
Universe’s infancy. Crucially, the attenuation of Ly$\alpha$ photons can probe
the presence of neutral hydrogen in the intergalactic medium (IGM) (e.g,
Miralda-Escudé, 1998; Malhotra & Rhoads, 2004; Dijkstra, 2014), but the
photons also undergo complicated resonant scattering within the galaxy,
complicating our understanding of how much of the emission exits the ISM and
circumgalactic medium (CGM) and enters the IGM in the first place. Recent
attempts to use Ly$\alpha$ as a reionization probe have struggled to account
for the intrinsic effects of host galaxy properties on the Ly$\alpha$
luminosity before the radiation encounters the IGM, leaving an unknown
systematic uncertainty present in their results. The most detailed
spectroscopic studies of post-reionization LAEs point to the covering fraction
of optically thick neutral hydrogen (e.g. Reddy et al. 2021) as the key
predictor of Ly$\alpha$ escape, but such observations remain expensive and
time intensive. Finding correlations between Ly$\alpha$ emission and global
properties such as mass and star formation activity, which photometry can
reliably measure even at very high redshifts, could be a path forward to
predicting galaxies’ intrinsic Ly$\alpha$ output.

Small LAE sample sizes ($<20$) were typical a decade ago, and although
recently large samples with $>$1000 objects have been amassed using narrow-
band surveys (e.g. Sobral et al. 2018, Ono et al. 2021), spectroscopically
confirmed samples remain small. This has statistically hindered the efficacy
of studies of global property correlations with Ly$\alpha$ emission. The
HETDEX project (Hill et al. 2008, Hill et al. 2021, Gebhardt et al. 2021) is
in the process of discovering a transformative sample of LAEs, clearing the
way for the community to obtain a better understanding of this intriguing
population. The un-targeted (targets not pre-selected), spectroscopically
selected HETDEX LAE sample at $1.9<z<3.5$ provides a unique vantage point of
galaxy evolution, as these galaxies probe the lower-mass end of the galaxy
distribution, making them analagous to typical galaxies discovered in the
epoch of reionization (e.g., Finkelstein et al. 2010).

As the first step toward realizing HETDEX’s ability to unlock LAEs as probes
of reionization, we present an initial study detailing how to link detections
from the survey to imaging counterparts, and we provide an SED fitting
analysis of their stellar population properties. Our modest sample of 72 LAEs
in the GOODS-N field will pave the way for future large samples from HETDEX to
obtain the best understanding of LAEs to date. In §2 we describe how we built
our sample and selected imaging counterparts. In §3 we describe our SED
fitting procedure. We present our results in §4, comparing them to other
studies, and we discuss our interpretations in §5. Finally, we attempt to
predict the Ly$\alpha$ emission from a sample of epoch of reionization (EoR)
galaxies in §6 and summarize this study in §7. In our analysis, we adopt a
flat $\Lambda$CDM cosmology with $H_{0}=70\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$ and
$\Omega_{\mathrm{m}}=0.30$.

##  2 Methodology

In order to explore how Ly$\alpha$ emission from galaxies depends on stellar
population properties, we built a sample of LAEs using emission line
detections from the HETDEX survey, carefully identifying them as Ly$\alpha$ or
other contaminant features, such as [O ii]$\lambda\lambda\ 3726,3729$, which
is unresolved at HETDEX resolution. We then created a procedure for assigning
the line detections to imaging counterparts in $HST$ data so that we could
proceed with fitting their SEDs.

###  2.1 The HETDEX Survey

With HETDEX the upgraded Hobby-Eberly Telescope (Ramsey et al. 1994, Hill et
al. 2021) is observing an area of $540\ \mathrm{deg}^{2}$ in the north
Galactic cap and on the celestial equator using up to 78 pairs of integral-
field spectrographs that span $350-550\ \mathrm{nm}$ at $R\sim 800$. Each
spectrograph pair is fed by an integral field unit (IFU) of 448
1.5$\arcsec$-diameter fibers which cover a 51$\arcsec$ $\times$ 51$\arcsec$
region on the sky with 1$/$3 fill factor (Kelz et al., 2014; Hill et al.,
2021). Each HETDEX observation consists of three 6-min dithered exposures to
fill in the area between fibers, each with $>$30,000 individual fibers. The
majority of these fibers just contain blank sky, but some subset contain
continuum sources such as stars or emission lines from both nearby and distant
galaxies.

Gebhardt et al. (2021) describe the data reduction and calibrations needed to
convert the raw observations into a three dimensional spectroscopic data set
as well as the methods used to detect emission lines contained in the millions
of observed spectra. As a brief summary, HETDEX reductions involve three types
of calibration frames: biases (taken nightly), pixel flats (taken yearly using
a laser-driven light source), and twilight sky flats (taken nightly and
averaged monthly), which are used for bias subtraction, bad pixel masking,
fiber profile tracing, wavelength calibration, scattered light removal,
spectral extraction, fiber normalization, spectral masking, and sky
subtraction. These frames, combined with sky background on science images,
produce a wavelength calibrated, sky-subtracted spectrum for each fiber in the
array.

Astrometric calibrations are achieved by measuring the centroid of each field
star and comparing their positions on the IFUs to the stars’ equatorial
coordinates in the Sloan Digital Sky Survey (SDSS; York et al., 2000;
Abazajian et al., 2009) and Gaia (Gaia Collaboration et al., 2018) catalogs.
This process typically results in global solutions which are good to $\sim
0.2\arcsec$ and no worse than $\sim 0.5\arcsec$, with the exact precision of a
measurement dependent upon the number of IFUs in operation at the time of the
observation.

To find emission lines, the data pipeline searched every spatial and spectral
resolution element in the internal HETDEX data release 2 (HDR2) to look for a
peak in signal. Regions of enhanced signal were fit with a single Gaussian
model with a constant continuum level, a model found adequate for potentially
asymmetric line profiles by Gebhardt et al. (2021) because of the low
resolution of the VIRUS spectrographs and low signal to noise ($S/N$) of
typical sources. The exact location was determined by rastering on a grid and
maximizing the line’s signal-to-noise. An internal catalog of high-quality
emission lines was generated by Mentuch Cooper et al. (in preparation), and we
drew our initial sample from the HDR2 version of that catalog. The catalog
reduced the raw detected line emission sources as described in Gebhardt et al.
(2021) into a more robust sample by passing the observations through a quality
assessment pipeline and limiting various fitted line parameters. Specifically,
emission lines were required to have a quality of fit, $\chi^{2}<1.2$ and a
linewidth, $\sigma$, in the Gaussian model between 1.7 Å and 8 Å. The full
HETDEX survey will eventually detect $\sim 1$ million LAEs, providing an
incredible opportunity to study such objects, but our analysis is focused on
LAEs discovered in 2018–2020 data from a HETDEX science verification field in
GOODS-N, a roughly $10^{\prime}\times 16^{\prime}$ field centered at (J2000)
$\mathrm{12^{h}36^{m}55^{s}},62^{\circ}14^{\mathrm{m}}15^{\mathrm{s}}$
(Giavalisco et al. 2004, Grogin et al. 2011, Koekemoer et al. 2011) because we
required deep, multi-band imaging to study each galaxy’s stellar populations.

###  2.2 Sample Selection

We visually inspected HETDEX detections in GOODS-N to obtain a clean sample of
LAEs. To get initial candidates, we applied various quality cuts to the
curated catalog for data release HDR 2.1.2 (Mentuch Cooper at al., in
preparation). We restricted emission line detections to those with signal-to-
noise ratio $S/N>5.5$ to limit the fraction of spurious detections from noise
fluctuations to less than 5% (see Gebhardt et al. 2021) as well as
$\chi^{2}<1.6$ for the Gaussian model fit, which was a value tuned to remove
the most obvious artifacts while retaining the largest sample for inspection.
We required emission line full-width at half-maximum (FWHM) between 3.4 Å and
24 Å, where the lower bound removed exceedingly narrow peaks arising from
unidentified cosmic rays and the upper bound removed emission generated by
broad-line AGNs, which we considered contaminants in this study (see §3.1). We
further only included observations with throughput $>0.07$ for reliable flux
measurements minimally affected by cloud cover, and seeing below 2.8′′ to
enable continuum counterpart identification. We did not remove “repeat”
detections coincident spatially and spectrally resulting from the survey
revisiting the field multiple times in order to ensure we found as many
Ly$\alpha$ detections as possible. We excluded data in GOODS-N taken prior to
2018 as they included significant artifacts from early CCDs that had been
replaced by 2018.

Finally, we did not initially remove any detections based on the Bayesian
probability values used to help determine the identity of an emission line as
Ly$\alpha$ vs [O ii]$\lambda\lambda\ 3726,3729$, such as P(Ly$\alpha$). These
probabilities, which are calculated by the HETDEX team based on the work of
Leung et al. (2017) and Farrow et al. (2021), leverage the inherent
differences between the emission line luminosity and equivalent width (EW),
$W_{\lambda}$, distribution functions of LAEs and [O ii] emitters to identify
single emission line detections using information about the line flux and
continuum emission, when available. During the process of visual inspection,
we used the statistic to guide our identifications, and we make
recommendations for using quality cuts based off this statistic at the end of
this section. We do not believe that keeping this statistic visible to the
classifier biased our results because we implemented an independent procedure
(see §2.3) to distinguish LAEs from low-redshift counterparts that relied on
SED fitting.

After applying quality cuts, we began with 842 detections (of which $\sim$500
were “unique” in the sense that there were no other emission line detections
within 3$\arcsec$ spatially and 6 Å spectrally). To inspect each detection, we
used the HETDEX Emission Line eXplorer tool ELiXer (Davis et al., in
preparation), which shows measured quantities for the emission lines such as
$S/N$, line width, line fit $\chi^{2}$, the continuum estimate, the Bayesian
probability for Ly$\alpha$ emission described above, as well as useful visual
information, such as cutouts of the 2D spectra for several fibers containing
the feature, the Gaussian model fit to the feature, the full 1D spectrum, and
any imaging and catalog data uploaded in the HETDEX pipeline.

We rated our confidence in a detection on a scale of 0-5 using a customized
widget tool that allows interactive classification of detected sources based
on the its Elixer Report (see Figure 1). Additionally, other classifications
include “artifact,” a false detection caused by a malfunction in the
instrument or the reduction pipeline, low-redshift sources, and “other” for
miscellaneous objects like meteors. To qualify for a classification of 4 or 5
(a high-confidence LAE by our definition), a detection had to meet the
following criteria:

  * •

A clear emission line in at least one fiber in the un-smoothed 2D spectrum, or
a probable emission line in at least two fibers. Since each point-spread-
function (PSF) covers multiple fibers (due to the dithering pattern), we
expected strong emission to be seen in more than one fiber, increasing the
likelihood of a real detection.

  * •

No obvious defects at the emission line location in the pixel flat or sky
subtraction cutouts. This eliminated hot pixels, sky model residuals, charge
traps, and other artifacts from the sample.

  * •

A Gaussian plus constant continuum model fit that adequately matched the data
and did not have a FWHM far below the spectral resolution of $\sim 6$ Å.

  * •

A line peak that exceeded the typical noise level in multiple pixels in the 1D
spectrum.

  * •

No source at the line’s detection position brighter than roughly $m_{AB}=24$
in the imaging cutouts, if available. The high equivalent widths of sources
fainter than this threshold drastically decrease the likelihood of
contamination by [O ii] emitters (see Figure 6 in Leung et al. 2017), though a
few low equivalent width, luminous LAEs can be missed with this requirement.

As the [O ii]$\lambda\lambda\ 3726,3729$ emission feature falls into the
$3500\ \mathrm{\AA}<\lambda<5500\ \mathrm{\AA}$ spectral range for $z<0.5$,
the imaging proved crucial in choosing between high-redshift LAEs and
interloping [O ii] emitting galaxies.

Figure 1: A section of an ELiXer report for detection ID 2100325857, the line
from an LAE in our sample. The report contains information about the detected
line as well as imaging at the detection position. The 2D spectra from the
four fibers contributing the highest $S/N$ to the detection are in the top
left corner; the stacked signal is shown on the top row, outlined in black.
The pixel flats and smoothed 2D spectra are displayed in the right two
columns. Plots to the right show the sky subtraction and the emission line
model fit (black line). The middle row contains $10\arcsec$ imaging cutouts
with the fiber positions shown in gray in the first image and the locations of
cataloged sources marked with colorful boxes in all subsequent images. White
boxes indicate sources too far to be considered, and white circles show the
aperture for brightness measurement in F606W. The AB magnitude of the nearest
likely source is reported. The top right corner shows the photo-$z$
probability distributions calculated by the CANDELS team (B. Andrews et al.,
in preparation) for sources of matching color in the imaging, and the [O ii]
and Ly$\alpha$ redshifts are shown as vertical dashed green and red lines.
Finally, the bottom row shows the full 1D spectrum. Flux densities have units
of $\mathrm{erg\ s^{-1}\ cm^{-2}\ 2\AA^{-1}}$

.

Figure 2: A section of an ELiXer report, in the same format as Figure 1, for
detection ID 2100037191, corresponding to the [O ii]$\lambda\lambda\
3726,3729$ feature in a galaxy at $z\approx 0.24$. The black trace in the 2D
fiber spectrum (blue rectangle) indicates a clear detection of continuum
emission, which is also evident in the 1D spectrum in at the bottom of the
figure. The imaging shows a large, bright source ($m_{AB}=22.0$) centered on
the detection position, and the source has a cataloged spectroscopic redshift
consistent with [O ii] indicated as an open blue circle in the photometric
redshift plot in the upper-right corner.

Figure 1 shows an example ELiXer report for a source classified as a high-
redshift LAE. Note that for readability, tabulated numeric information such as
P(Ly$\alpha$), line flux, line model $\chi^{2}$, and more was cropped out of
this visualization, but was visible to the classifier. In Figure 1, A clear
emission feature is present as a black signal in three out of the four 2D un-
smoothed fiber spectra, the sky subtraction looks clean, the model fit
accurately represents the data, and the image stamps show a number of faint
sources with photometric redshift estimates reasonably close to the Ly$\alpha$
redshift (shown by the vertical red dashed line). Figure 2 shows a clear
example of a low-redshift object detected by its [O ii] emission line. As in
Figure 1, the line appears strong in multiple fibers, and the sky subtraction
and model fit present no concerns. Characteristically of a brighter low-
redshift galaxy, continuum emission is visible as a horizontal black trace in
the fiber spectra, and a large, bright object appears in the $HST$ image
stamps. In this case, the object is in fact a cataloged [O ii] emitter, but
even without such information this would be a clear low-redshift
classification. In both of these cases, no other emission lines are detected,
or would be expected to be detectable, across the observed wavelength range.

After classifying each detection, we obtained $\sim 200$ detections
categorized as high confidence LAEs (scores of 4-5) and almost three times as
many classified as low-$z$ sources (Figure 3). Note that we did not include
detections with scores of 3 or below for initial study as we want the cleanest
sample possible. To assess the HETDEX collaboration’s built-in Bayesian
classification probability, P(Ly$\alpha$), we plotted that statistic for all
of our detections classified as either low-$z$ galaxies or LAEs. Figure 4
shows that true LAE detections rarely score low in the P(Ly$\alpha$)
statistic, but a few low-$z$ sources can score in the intermediate range. For
this reason, we suggest future studies can dramatically reduce the amount of
visual inspections needed by adopting a cutoff of P(Ly$\alpha$)$\gtrsim 0.6$
for LAE candidates.

Figure 3: The distribution of visual classifications of candidate detections.
Making cuts based on P(Ly$\alpha$) to remove low-z sources can dramatically
reduce the visual inspection workload. Figure 4: The distribution of
P(Ly$\alpha$) for detections visually classified as low-redshift galaxies
(light blue) and high-confidence (scores of 4 and 5) LAEs (light red).
Adopting a minimum threshold for P(Ly$\alpha$) can remove a large fraction of
low-redshift interlopers without eliminating very many LAEs.

To finalize our sample, we removed detections of the same source (since the
GOODS-N field was observed multiple times between 2018–2020), by selecting the
highest $S/N$ measurement of all detections grouped within 2$\arcsec$ and one
spectral resolution element (6 Å). Our final emission line sample consisted of
94 high-confidence Ly$\alpha$ detections (with classification scores of 4-5).

###  2.3 Counterpart Identification

In order to study the stellar populations of the LAEs in our sample, we
developed a method to match the un-targeted spectroscopic detections to
counterparts in $HST$ imaging of the GOODS-N field.

The overall astrometric precision of a HETDEX observation is $\sim$0.2″.
However, due to the 1.5″ diameters of the fibers, the typical seeing, and the
3-dither pattern, the position of an individual (faint) LAE is known to no
better than $\sim$0.5″. Since the HST images have a resolution that is
$\sim$20 times higher than this, great care is needed to ensure an emission-
line source is matched with the correct counterpart.

We used the imaging obtained by the Great Observatories Origins Deep Survey
(GOODS Giavalisco et al., 2004) with the optical ACS camera, and the Cosmic
Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS Grogin et
al., 2011; Koekemoer et al., 2011) with the WFC3/IR infrared camera, using the
internal CANDELS team’s reduced mosaics for each filter. This dataset consists
of imaging in nine filters (F435W, F606W, F775W, F814W, F850LP with ACS, and
F105W, F125W, F140W and F160W with WFC3/IR). We made use of the photometric
catalogs derived by Finkelstein et al. (2021) which used SourceExtractor
(Bertin & Arnouts, 1996) in two-image mode to create a F160W-selected catalog,
coupled with using the Tractor (Lang et al., 2016) to perform deblended
photometry on the deep S-CANDELS (Ashby et al., 2015) Spitzer/IRAC 3.6 and 4.5
$\mu$m imaging. Further details on the cataloguing process are available in
Finkelstein et al. (2021). Similar to the widget used to classify detections
as Ly$\alpha$ , we created a visual inspection tool that provided information
about the distance between the centroid of the HETDEX emission location and a
given imaging source, the HETDEX emission-line strength when re-extracted
centered at the imaging counterpart position, and the goodness of an SED fit
assuming the Ly$\alpha$ redshift, $z_{\mathrm{Ly\alpha}}$.

Before selecting counterpart candidates, we optimized our search by developing
a deep photometric catalog using a stacked image across all $HST$ filters in
GOODS-N. Each pixel value in this image and its error was computed using an
inverse variance weighted average across $N=9$ filters with pixel value
$p_{i}$ and rms error $\sigma_{i}$ given by Equation 1.

$\bar{p}=\frac{\sum_{i}^{N}p_{i}\sigma_{i}^{-2}}{\sum_{i}^{N}\sigma_{i}^{-2}}\
,\ \sigma_{\bar{p}}=\left(\sum_{i}^{N}\sigma_{i}^{-2}\right)^{-1/2}$ (1)

Since LAEs are often low-mass, faint systems, this stacked image improved our
chances of identifying the continuum source corresponding to the detected
emission line.

We then used SourceExtractor (Bertin & Arnouts, 1996) to detect the faintest
possible sources in the stacked image, requiring a source to have 5 contiguous
pixels with $S/N$$>1.6$. Following the procedures outlined in Finkelstein et
al. (2021) we used the same software in two-image mode to measure the flux in
each filter and applied the appropriate aperture correction obtained from
simulations. We performed extinction corrections using a Cardelli extinction
law with $R_{V}=3.1$ for the Milky Way (Cardelli et al., 1989). We then
compared the fluxes measured in this catalog to the F160W-selected catalog of
Finkelstein et al. (2021) and found the flux measurements to have no
systematic offset and minimal scatter. Figure 5 shows the fractional error of
the stacked catalog photometry compared to the Finkelstein et al. 2021
photometry as a function of source brightness in the $I$-band. The median
offset is zero with scatter of roughly 25% for fluxes near 100 nJy, in
agreement with the typical error bars for such sources, providing confidence
in the fidelity of the stacked catalog. In all subsequent analysis, we
defaulted to using measurements from the Finkelstein et al. (2021) catalog for
sources detected in both, and we only used photometry from the stacked catalog
for five LAEs in our sample unique to it.

Figure 5: The fractional error between $I$-band (F775W) fluxes of sources
measured in our derived stacked-detection catalog matched to sources in the
original catalog from Finkelstein et al. (2021) as a function of flux in the
latter catalog. The median offset is indicated with the dashed pink line, and
its value is given with text, showing good agreement between these two
catalogs.

After generating the catalog from stacked imaging, we identified all imaging
sources within $3\arcsec$ of the HETDEX detection position as possible LAE
counterparts. Since the typical image quality of the HETDEX observations used
here has a point spread function (PSF) of $\sim 1.7\arcsec$, the $3\arcsec$
annulus served as a generous aperture around the Ly$\alpha$ centroid to
encompass all possible counterparts for the detected emission.

Figure 6: An example of the visualizations used to select imaging counterparts
for Ly$\alpha$ detections. The figure shows a 10$\arcsec$ image cutout with
the detection positioned marked with a magenta star and sources marked with
“X”s of various colors, 1D spectra extracted at each source position, the
original detection 1D spectrum, and SED fits (with redshift fixed assuming
Ly$\alpha$ emission) of all sources with significant measured fluxes. The
table in the top-right contains the $\chi^{2}$ value of each SED fit, the
separation between each source and the emission line detection position
(labeled $\theta$), and the significance (labeled SNR) and the line flux in
$\mathrm{erg\ s^{-1}\ cm^{-2}}$ of the emission line extracted at the source
position. The colors are consistent across all plots and tables, so each
source corresponds to a unique color. The SED plots also contain normalized
filter response curves as gray lines. In this case, while the red and orange
sources have similarly small distances from the detection position and similar
line fluxes, the SED fit $\chi^{2}$ strongly favors the orange source to be an
LAE at $z=2.90$.

We selected imaging counterparts based on the neighboring sources’ angular
distances from the detection, significance of emission extracted at the source
positions, and goodness of SED fits performed by fixing the redshift assuming
a Ly$\alpha$ detection. First, we measured the on-sky angular separation from
the detection position to the position of each possible source in the
photometric catalog (labeled $\theta$ in Figure 6). Then, for each source, we
used the HETDEX API script,
get_spectrum.py111https://github.com/HETDEX/hetdex_api, to perform an
aperture-weighted optimized spectral extraction (following Horne 1986) at the
source position to obtain a 1D spectrum. We created a Markov chain Monte Carlo
(MCMC) line-fitting code using emcee (Foreman-Mackey et al., 2013) to fit a
model to the feature to estimate its flux and significance. Our model
consisted of two components: a linear trend with slope $m$ and intercept $b$,
which captured any underlying continuum, and a Gaussian with total flux $F$
and standard deviation $\sigma$ to fit the line profile.

$f_{\lambda}=m(\lambda-\lambda_{0})+b+\frac{F}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(\lambda-\lambda_{0})^{2}}{2\sigma^{2}}\right]$
(2)

In the model, $\lambda_{0}$, the wavelength of the emission line, was allowed
to vary by $\pm$ one pixel (2 Å) from the detection wavelength reported by
HETDEX. For each fit, we measured an effective $S/N$ ratio (labeled SNR in
Figure 6) by comparing the median value of the line flux to the standard
deviation of the line flux for the last 20% of the MCMC sampling chain, which
had converged at that stage of sampling. To limit computation time for future
counterpart identification steps, we ruled out any counterpart candidates that
did not have any indication ($S/N>1$) of an emission feature at the pixel
corresponding to the detected wavelength. Finally, for those sources with
significant emission, we performed SED fitting with Bagpipes (see Section 3.1
for a full description of this procedure), fixing the redshift as $z_{{\rm
Ly}\alpha}$. Our simple SED model for counterpart identification included free
parameters for stellar mass, metallicity, dust extinction, and SFH, and we
adopted the Calzetti et al. (1994) dust attenuation law, the Chabrier (2003)
initial mass function, and a delayed-$\tau$ SFH. At this stage, we did not
include any IRAC fluxes in our fits since those fluxes depend sensitively on
deblending, which is unreliable when sources are crowded. Furthermore, between
$1.9<z<3.5$, there are no strong spectral features at the rest-frame
wavelengths probed by IRAC, and redshift-sensitive features such as the 4000 Å
break are adequately covered by $HST$. We then visually inspected the
separations, spectral extractions, and SED fits of all candidate counterparts
to choose the one most likely to be the detected LAE.

Figure 6 shows an example of our approach. Separate sources are marked with an
‘‘X’’ and the color of the mark corresponds to the color of the table row,
spectrum, and SED in the subsequent plots. In this case, the red and and
orange sources within $0\farcs 5$ of the detection position (magenta star)
show similar extracted emission line flux at the detection wavelength.
Crucially, the SED fit for the red object poorly matches the data when fixing
the redshift as $z_{Ly\alpha}$, but the orange object has a fit in excellent
agreement with its observed SED based on the $\chi^{2}$ statistic. Therefore,
in this example case we selected the orange object as the detected LAE at
$z=2.90$. We followed the same process to identify counterparts for the other
Ly$\alpha$ lines in our sample.

By studying the distribution of our counterparts in the parameter space of
separation, signal-to-noise of emission, and SED $\chi^{2}$, we found no
obvious way to select counterparts reliably based on these numbers alone, but
we did find favorable regions. Figure 7a shows the distribution of separation
from the detection positions for sources we identified as LAEs and sources
that just happened to be nearby. Clearly, it was exceedingly unlikely that the
true counterpart lay farther than 1$\arcsec$ away on-sky. For this reason, we
could very reasonably shrink our selection criteria from all sources within
3$\arcsec$ to roughly 1$\arcsec$ without significant loss of LAEs. In terms of
emission line $S/N$ (compared to the measured value of the detection itself),
we found that, while typically the identified counterparts had stronger
emission, the HETDEX PSF caused the extracted flux to not depend sensitively
enough on position to clearly identify the counterpart for sources separated
by less than 1$\arcsec$. This is clearest in Figure 8a, which shows that true
counterparts and close neighbors show overlap in the $S/N$, separation plane.
Note that the different on-sky centroids for emission line extraction between
the counterparts and the original detection allow for the values of the $S/N$
ratios in Figure 7b to be greater than unity. Finally, we note that, while
most of the LAEs in our sample had $\chi^{2}$ values in good agreement with
the $z_{\mathrm{Ly\alpha}}$ hypothesis, many neighboring galaxies also had low
$\chi^{2}$, as shown in Figure 7c. We attribute the low $\chi^{2}$ values for
non-counterparts to our inclusion of such faint objects, which have large flux
errors and are thus easily fit by a wide range of models.

Figure 7: The distributions, expressed as the fraction of objects in a given
bin, of (a) separation, (b) position-extracted emission line signal-to-noise
(relative to that of the detection), and (c) SED $\chi^{2}$ assuming
$z_{\mathrm{Ly\alpha}}$, for objects identified as the detection imaging
counterparts and those that happened to be spatially coincident. Histograms
are normalized to the population size. The top panel indicated that finding a
counterpart with an imaging separation larger than 1$\arcsec$ from the
detection position is exceedingly rare. Figure 8: (a) The 2D distribution (in
$S/N$ and separation space) for objects identified as the imaging counterparts
for emission line detections and those that happened to be spatially
coincident. (b) The same plot in SED $\chi^{2}$ and separation space. These
two figures show substantial overlap in these parameter spaces for true LAEs
and neighboring sources, motivating the benefits of detailed visual inspection
shown in Figure 6.

After visually vetting all detections in our sample of 94 Ly$\alpha$ lines, we
found 6 instances of detected emission with no continuum-detected counterpart.
Since we could not study the properties of an LAE without photometry, we
removed these objects from the final analysis. Furthermore, we removed 16
objects from the sample due to the following quality concerns. We eliminated
the LAE corresponding to HETDEX detection ID 2100245124
(RA,DEC=189.346621∘,62.260662∘) from our sample as it was the only counterpart
with an X-ray detection in the catalog of Xue et al. (2016), indicating the
galaxy hosted an AGN. Since our SED fitting code did not have an AGN template,
we could not reliably report the physical properties of this object. We also
eliminated the detection for ID 2100171783, as the counterpart inspection
revealed the Ly$\alpha$ emission line came from two probable LAEs separated by
less than $0\farcs 5$ meaning we could not assign flux accurately to each
source. Finally, we only analyzed objects detected in the $H$-band (F160W) of
the $HST$ imaging as well as at least two bluer bands in order to span the
rest-frame 4000 Å break at the sample redshift range. This serves as a crucial
feature for constraining galaxy masses and ages with SED fitting (e.g., see
Shapley et al. 2003). These choices limited our final sample size to 72 LAEs
in GOODS-N spanning $1.98<z<3.48$. For 5 of these objects, photometry was not
present in the catalog of Finkelstein et al. (2021), so we used photometry
from the stacked catalog described in §2.3. Appendix B, Figure 21 shows all
the HETDEX emission lines for the LAEs in the final sample, and Figure 22
shows $HST$ imaging in F160W ($H$-band) for all objects.

##  3 Analysis

After connecting HETDEX emission line detections with $HST$ imaging
counterparts, we leveraged SED fitting to measure the galaxies’ stellar
population properties. From the SED fits and emission line detections, we also
inferred the UV-slope and Ly$\alpha$ equivalent width.

###  3.1 SED Fitting with Bagpipes

We fit all LAEs in our final sample with Bagpipes (Carnall et al., 2018), a
flexible python code that rapidly generates galaxy model spectra through
stellar population synthesis using the 2016 version of the Bruzual & Charlot
(2003) stellar spectral libraries. It explores the high-dimensional, multi-
modal, and degenerate (e.g., age-dust-metallicity) model parameter space using
the MultiNest algorithm (Feroz & Skilling, 2013).

Our sample in GOODS-N had photometry across nine $HST$ filters ranging from
0.4 to 1.6 $\mathrm{\mu m}$ as well as two $Spitzer$/IRAC channels centered at
3.6 $\mathrm{\mu m}$ and 4.5 $\mathrm{\mu m}$. Translating to the rest-frames
of the objects in our sample at $1.9<z<3.5$, these filters probed the UV,
optical, and near-infrared (NIR) energy output of our objects.

The filter coverage of our sample of LAEs motivated our choice of SED modeling
parameters. Table 1 shows the names and units of the free parameters in our
model, as well as the prior probability distributions assumed in our Bayesian
framework. We adopted a delayed-$\tau$ SFH, defined as:

$\mathrm{SFR}(t)\propto\begin{cases}(t-t_{0})e^{-(t-t_{0})/\tau}&t>t_{0}\\\
0&t<t_{0}\end{cases}$ (3)

This flexible SFH allows for star formation to be either rising, peaking, or
falling, as opposed to the common exponentially declining model that only
allows for falling SFRs over time. For example, Lee et al. (2010) found that
SED fitting that adopted rising SFHs matched the stellar masses and SFRs from
semi-analytic models for galaxies at $3<z<6$ better than exponentially
declining models, while Papovich et al. (2011) found similar results favoring
rising SFHs for real galaxies at $z=$ 4--7. We fit the $e$-folding scale of
the SFH, $\tau$, the age of the Universe at the onset of star formation,
$t_{0}$, the stellar mass formed, $M_{\mathrm{form}}$, the global metallicity,
$Z$, the dust extinction in the $V$-band, $A_{V}$, and the ionization
parameter, $\log U$, defined as the log of the ratio of the number densities
of ionizing photons and hydrogen atoms. Though we fit the total stellar mass
formed by a galaxy, $M_{\mathrm{form}}$, we report its stellar mass at the
redshift of observation excluding remnants, and we denote that stellar mass
$M_{\star}$. We note that some of the parameters (namely $Z$ and $U$) are not
expected to be well-constrained by our photometric data. Nonetheless we allow
them to vary within our imposed priors such that the uncertainties in the
other parameters include the uncertainties in these parameters. We adopted the
Calzetti et al. (1994) dust attenuation law for star-forming galaxies and the
Chabrier (2003) initial mass function.

Parameter | Prior | Bounds | Units
---|---|---|---
$t_{0}$ | Uniform |  0, $T(z)$ | Gyr
$\tau$ | Uniform | 0.3, 10. | Gyr
$M_{\mathrm{form}}$ | Log Uniform |  $10^{6}$, $10^{12}$ | $\mathrm{M_{\odot}}$
$Z$ | Log Uniform |  $10^{-5}$, 2 | $\mathrm{Z_{\odot}}$
$A_{V}$ | Uniform | 0, 2 | mag
$\log U$ | Uniform | -4, -2 | -
Table 1: Free parameters and their prior probability distributions for SED
fitting. In our galaxy models, the redshift, $z$, was fixed based on the
observed wavelength of Ly$\alpha$ from HETDEX. $T(z)$ refers to the age of the
Universe at redshift $z$. Note that we fit the cumulative stellar mass formed,
$M_{\mathrm{form}}$, from which the stellar mass (excluding remnants) at the
object redshift was computed within the Bagpipes (Carnall et al., 2018) code.

All 11 filters were not necessarily included for every galaxy SED fit in our
sample. For example, due to the large PSF of the IRAC imager, modeling sources
in crowded fields of view and deblending the flux contribution of each source
is crucial to accurately measuring the NIR fluxes of our LAEs. Although the
catalog we used performed deblended photometric modeling with the IRAC PSF,
this process can fail in crowded regions. We thus visually inspected all IRAC
residual maps for objects in our sample and removed the IRAC fluxes from our
SED fitting if there were obvious problems in the deblending procedure. For
the 5 objects not present in the catalog of Finkelstein et al. (2021), we did
not have IRAC measurements. Furthermore, because the purpose of our analysis
was to study the SED-derived properties of our LAEs in relation to their
Ly$\alpha$ emission, we did not want Bagpipes’s modeling of Ly$\alpha$
emission or the IGM attenuation to bias our results. For this reason, we
masked out all filters whose bandpass extended blue-ward of the observed
Ly$\alpha$ line; thus the $B$-band (F435W) and sometimes the $V$-band (F606W)
was excluded, depending on redshift.

Figure 9: (Left): An example fit to the SED of an LAE in the sample. Data are
shown by blue circles, and the 68% spread for the posterior model photometry
and spectrum are shown by the green rectangles and orange shaded lines,
respectively. A maroon open circle indicates the measured flux for a filter
masked during SED fitting. In this case, the $B$-band was masked out since it
includes the Ly$\alpha$ emission line. Our imaging data constrain the $4000\
\mathrm{\AA}$ break and rest-optical colors. (Top-right): $5\arcsec$ square
image cutout for the source in the F160W $HST$ filter. The pink cross
indicates the source position, and the blue plus sign and dashed blue circle
indicate the detection position and the FWHM of the HETDEX fiber PSF. (Bottom-
right): 1D extracted spectrum for this source, centered on an $80\
\mathrm{\AA}$ window around the Ly$\alpha$ emission line. The solid green line
indicates the HETDEX Gaussian model fit to the data.  Figure 10: A ‘‘corner’’
plot of the fit in Figure 9 for object ID 10388. The 1D histograms are shown
on the diagonal for the posterior distribution of each free parameter in our
model (see Table 1). The 2D histograms show the correlations of all parameters
with one another, where contour lines are drawn for each $\sigma$ level. With
our broadband photometry, we constrained ages, masses, and dust extinctions
well.

Figure 9 shows an example Bagpipes SED fit for an LAE in our sample. We
plotted the $1\sigma$ spread on the model photometry as rectangles as well as
the $1\sigma$ spread on the underlying model spectrum computed by evaluating
the 16th and 84th percentiles of the posterior models. In this example, the
fit did an excellent job matching salient features like the rest-frame 4000 Å
break and nebular emission in the rest-frame optical region. We estimated
galaxy properties using the posterior distributions for all free parameters
explored by Bagpipes. Figure 10 shows an example ‘‘corner’’ plot (produced via
Foreman-Mackey 2016), where all free parameters are plotted against each other
for easy assessment of constraints and correlations. Stellar mass, time since
the onset of star formation, and dust extinction were constrained well, while
metallicity, $\tau$, and ionization parameter were not well-constrained by our
broadband photometry data. Figure 20 in Appendix B shows the SED fits for all
72 LAEs in our final sample.

###  3.2 Measuring $W_{\lambda}({\mathrm{Ly\alpha}})$ and $\beta$

Emission line strengths can be represented by the parameter equivalent width
(EW or $W_{\lambda}$), which represents the width of a rectangle drawn to the
same height as the continuum needed for the rectangular area to match the area
under the emission line. To estimate the equivalent width of HETDEX Ly$\alpha$
detections, we used the measured line flux and error from the internal HDR
2.1.3 catalog computed by optimally extracting flux from all fibers within a
3.5 ″ radius circular aperture (roughly 15-20 individual fiber spectra)
contributing to the emission line detection (following Horne 1986), weighted
by the PSF of a point-source. We approximated the continuum flux density using
the Bagpipes sampled model spectra from the SED fit. We took the continuum
flux density to be the median value of all 500 sampled spectra averaged
between 1250 and 1300 Å in the given object’s rest-frame, and we computed the
$1\sigma$ error using half the spread between the 16th and 84th percentiles of
those values. This method allowed us to take advantage of complex computations
performed by Bagpipes to get a statistically representative estimate of the
continuum flux density instead of using a coarse approximation based off the
flux in one of our photometric bands. We evaluated the Ly$\alpha$ flux and the
continuum flux density in the observer-frame and translated to the galaxy
rest-frame by dividing by a factor of $(1+z)$ using the detected wavelength of
Ly$\alpha$.

$W_{\lambda}({\mathrm{Ly\alpha}})=\frac{F_{\mathrm{Ly\alpha}}}{f_{\lambda}}(1+z)^{-1}$
(4)

We measured $\beta$, the UV continuum slope (un-corrected for dust), using the
model spectra for galaxies in our sample following the method described in
Finkelstein et al. (2012). We masked the stellar and interstellar absorption
features in the rest-frame UV using the windows provided by Calzetti et al.
(1994), and we fit a linear model to the spectrum in log space ($\log
f_{\lambda}=\beta\log\lambda+C$) using polyfit from the python package Numpy
(Harris et al., 2020). We determined $1\sigma$ uncertainties on $\beta$ for
each object by measuring the distribution of values fitted to 500 spectral
models sampled from the posterior by Bagpipes.

##  4 Results

We measured various physical properties of objects in our LAE sample using the
posterior distributions returned by Bagpipes’ exploration of the parameter
space. We took the 16th and 84th percentiles of the posterior distributions to
represent the error bars on physical properties. Examples of such measurements
are shown in Figure 10 for a representative LAE in our sample.

###  4.1 SED-Derived Properties

Figure 11 shows the 1D distributions of posterior median values of stellar
mass ($M_{\star}$), star formation rate (SFR), specific star formation rate
(sSFR), dust extinction ($A_{V}$), mass-weighted age, and UV-slope ($\beta$)
for all objects in our final LAE sample. We found the median stellar mass of
our HETDEX LAEs to be $0.8^{+2.9}_{-0.5}\times 10^{9}\ \mathrm{M_{\odot}}$.
This stellar mass value lies near the median masses of LAEs selected in
narrowband imaging surveys covering redshifts similar to this study (e.g.
Guaita et al. 2011, Gawiser et al. 2007, Vargas et al. 2014, Kusakabe et al.
2018, Santos et al. 2020) and well below typical masses of Lyman-break
selected objects (e.g. Shapley et al. 2003, Papovich et al. 2001, Trainor et
al. 2019), which often have minimum masses an order of magnitude larger due to
the depth of the broadband imaging used in their selection.

We used our SED fitting procedure to obtain the attenuation in the $V$-band of
starlight due to dust for galaxies in the sample, and we obtained a median
value of $A_{V}=0.3^{+0.4}_{-0.1}$ mag. The presence of dust has been measured
in many other samples of LAEs with values of $A_{V}$ or $E(B-V)$ often falling
within a factor of two of this study (e.g. Guaita et al. 2011, Finkelstein et
al. 2009, Hathi et al. 2016, Kusakabe et al. 2018, Matthee et al. 2021).

Figure 11: The distributions of posterior median values for (clockwise from
top left) stellar mass, SFR, sSFR, dust extinction (in $V$-band mag), mass-
weighted age, and UV-slope for all objects in the sample. The 16th, 50th, and
84th percentiles are indicated by vertical dashed grey lines, and their values
are indicated with text in the same units as the x-axis labels. The LAEs in
our sample exhibit average properties similar to other LAE samples compiled at
comparable redshifts using narrow band selection.

Similar to dust reddening, our LAE sample has similar ages and star formation
rates to LAE samples in the literature compiled using narrow band or continuum
selection methods. Our SED-derived mass-weighted ages, typically spanning
0.05-0.5 Gyr, broadly agree with the narrow band samples of Acquaviva et al.
(2011), Finkelstein et al. (2009), Gawiser et al. (2007), and Vargas et al.
(2014). Our median SFR, $4.8^{+10.4}_{-3.8}\mathrm{M_{\odot}/yr}$, falls near
values reported by Gawiser et al. (2007), Hathi et al. (2016), and Kusakabe et
al. (2018), but falls 1 dex above the median SFR for LAEs found in the MUSE
HUDF Survey (Feltre et al., 2020). This discrepancy does not surprise us since
the MUSE HUDF LAE sample had a median mass roughly 0.5 dex lower than this
study, and their sample spanned $2.9<z<4.6$, probing an era of lower star
formation activity in the Universe than the one studied here (see Madau &
Dickinson 2014).

Our model included stellar metallicity and ionization parameter as free
parameters, but our broadband photometric data could not constrain those
values precisely (see Figure 10), since reliable estimates typically require
sensitive emission line diagnostics (e.g. Reddy et al. 2021), which were
coarsely probed at best by our filter set. For this reason, we do not present
or discuss our galaxies’ metallicities or ISM ionization conditions, but we
note that by letting these parameters vary, our posterior constraints on all
other parameters include the uncertainties in these quantities.

###  4.2 $W_{\lambda}({\mathrm{Ly\alpha}})$ Distribution

The equivalent width distribution of LAEs has been modeled by various authors
as exponential with the form given by Equation 5 (e.g. Gronwall et al. 2007,
Guaita et al. 2010, Wold et al. 2014, Jung et al. 2018).

$\frac{\mathrm{d}N}{\mathrm{d}W_{\lambda}}\propto e^{-W_{\lambda}/W_{0}}$ (5)

We show our sample’s rest-frame $W_{\lambda}({\mathrm{Ly\alpha}})$
distribution in Figure 12 with an $e$-folding scale $W_{0}=100\ \mathrm{\AA}$
drawn for comparison. We cannot measure the underlying distribution for LAEs
from our sample since we have not measured the completeness as a function of
equivalent width (which is complex due to our method of sample creation, and
not crucial for our study of stellar population properties). Various other
studies have precisely measured the Lyman-alpha equivalent width distribution,
such as Gronwall et al. (2007), who found an $e$-folding scale of
$76^{+11}_{-8}\ \mathrm{\AA}$ for a deep, narrow-band (NB) selected LAE sample
at $z=3.1$, Guaita et al. (2010), who measured $W_{0}=50\pm 7\ \mathrm{\AA}$
for a NB sample at $z=2.1$, and recently Santos et al. (2020), who measured
$W_{0}=129\pm 11\ \mathrm{\AA}$ for the full SC4K sample at $2<z<6$. We plot
some of these measured distributions in Figure 12 for comparison. It is
apparent that our sample becomes increasingly incomplete at EW $\lesssim$ 50
Å, due to a combination of the HETDEX flux limit, the emission-line
identification process, and our counterpart selection process.

Figure 12: The equivalent width distribution of LAEs in the sample. An
exponential distribution with $W_{0}=100\ \mathrm{\AA}$ is drawn in red for
comparison, as well as models fit by Gronwall et al. (2007) and Santos et al.
(2020). Our data favor models with larger values of $W_{0}$ to best match the
number of high-EW sources.

###  4.3 Correlations between $W_{\lambda}({\mathrm{Ly\alpha}})$ and Galaxy
Properties

We combined our SED-derived galaxy properties with the
$W_{\lambda}({\mathrm{Ly\alpha}})$ measurements described above in order to
assess correlations between Ly$\alpha$ emission and global galaxy properties.
We used $W_{\lambda}({\mathrm{Ly\alpha}})$ as a proxy for the fraction of
photons emitted as Ly$\alpha$ as opposed to $L_{\mathrm{Ly\alpha}}$, for
example, because the equivalent width more closely probes the physics
governing Ly$\alpha$ escape, whereas the flux also includes physics related to
the Ly$\alpha$ production rate. Figure 13 shows $M_{\star}$, specific star
formation rate (sSFR), star formation rate (SFR), dust extinction ($A_{V}$),
mass-weighted stellar population age, and UV-slope ($\beta$) plotted against
each galaxy’s $W_{\lambda}({\mathrm{Ly\alpha}})$ measurement. In the figure,
error bars denote the 16th to 84th percentile range, and we indicate Pearson’s
linear correlation coefficient, $r_{p}$, and its significance ($p-$value) with
text.

Figure 13: The relationship between $W_{\lambda}({\mathrm{Ly\alpha}})$ and
(clockwise from top left) stellar mass, SFR, sSFR, dust extinction (AV), mass-
weighted age, and UV-slope ($\beta$) for all objects in the sample. Grey lines
indicate the $\pm 1\sigma$ error bars on physical properties and equivalent
width. Open red circles show medians and standard deviations for properties in
equivalent width bins having equal numbers of objects. Pearson correlation
coefficients, $r_{p}$, and $p$-values are indicated for each plot. Stellar
mass and SFR exhibit the strongest correlations with
$W_{\lambda}({\mathrm{Ly\alpha}})$, while age and sSFR correlate moderately.
Surprisingly, no strong correlation exists with dust extinction.

Stellar mass and star formation rate both correlate strongly with
$W_{\lambda}({\mathrm{Ly\alpha}})$, with low mass, low SFR systems achieving
larger $W_{\lambda}({\mathrm{Ly\alpha}})$ than higher mass systems. The
correlation with mass has been established in the literature from studies of a
wide variety of galaxies such as LBGs, oELGs, and LAEs. It was noticed early
by Ando et al. 2006 and measured recently by many works such as Du et al.
2018, Marchi et al. 2019, Oyarzún et al. 2017, and Shimakawa et al. 2017.
Specifically, Weiss et al. (2021) found a negative correlation between the
Ly$\alpha$ escape fraction, $f_{esc}^{\mathrm{Ly\alpha}}$, and stellar mass
using data from the HETDEX survey. Additionally, Khostovan et al. (2021) found
an intrinsic, negative correlation between H$\alpha$ equivalent width and
galaxy stellar mass from a NB survey at $z\sim 5$. While the lack of low-EW,
low-mass systems can be driven by selection incompleteness, we should be
complete to high-mass, high-EW systems, yet these are seemingly rare.

Notably, our results show no significant anti-correlation between
$W_{\lambda}({\mathrm{Ly\alpha}})$ and dust extinction ($A_{V}$), whereas
numerous other studies of Lyman-alpha emission measured a clear relationship
that indicates dust hinders the ability of the Ly$\alpha$ photon to escape the
galaxy. For example, Shapley et al. (2003), Guaita et al. (2011), Du et al.
(2018), Hathi et al. (2016), Huang et al. (2021), Marchi et al. (2019),
Matthee et al. (2016), Reddy et al. (2021), Trainor et al. (2019), and Weiss
et al. (2021), all showed that dustier galaxies exhibit weaker Lyman-alpha
emission measured as $W_{\lambda}({\mathrm{Ly\alpha}})$ or have smaller
$f_{esc}^{\mathrm{Ly\alpha}}$. However, the lack of a significant anti-
correlation may be due to our limited sample size and small dynamic range in
dust attenuation. Moreover, objects with significant amounts of dust that
suppress their Ly$\alpha$ fluxes would not become members of our science
sample in the first place. The majority of our sample has $A_{V}<$ 0.3. We do
observe multiple galaxies with $A_{V}>$ 0.5, and interestingly these do not
all have low $W_{\lambda}({\mathrm{Ly\alpha}})$, implying that Ly$\alpha$ can
escape even from modestly dusty galaxies, which could indicate enhanced escape
due to outflows (e.g., Steidel et al. 2010, Erb et al. 2012) or a multi-phase
ISM (e.g. Finkelstein et al. 2009, Neufeld 1991).

Our Pearson correlation coefficient suggests a moderate correlation between
$W_{\lambda}({\mathrm{Ly\alpha}})$ and galaxy stellar mass-weighted age
($r_{p}=0.32$) in the sense that older galaxies exhibit larger
$W_{\lambda}({\mathrm{Ly\alpha}})$. Marchi et al. (2019) found a similar
result, obtaining a Spearman rank correlation coefficient of 0.40. This
contrasts with Pentericci et al. (2009) and Pentericci et al. (2010) who found
no strong dependence of Ly$\alpha$ equivalent width on age for LAEs and LBGs,
as well as Reddy et al. (2021) who found a weak negative correlation between
the two measurements for star-forming galaxies in the same redshift range
probed by this study.

Finally, a moderate negative correlation exists between sSFR and
$W_{\lambda}({\mathrm{Ly\alpha}})$, though the large error bars for our
measurements of sSFR weaken the reliability of the correlation. For
comparison, Hathi et al. (2016) found no significant correlation between the
two properties for a sample including Ly$\alpha$ in absorption and emission.

We also plot SFR against $M_{\star}$ for all objects in our sample in Figure
14 to see how our galaxies compare to other objects at similar redshift in
relation to the star-forming main sequence (SFMS). We include the best-fit
line found by Sanders et al. (2018) for star forming galaxies in the MOSDEF
survey at $z\sim 2.3$. Note that masses derived for that study used the
Chabrier (2003) IMF and Calzetti et al. (2000) dust curve but stellar
population synthesis models from Conroy et al. (2009). We also use a colorbar
to show the value of $W_{\lambda}({\mathrm{Ly\alpha}})$ for each galaxy. The
position of LAEs on the SFMS remains somewhat controversial. Studies such as
Vargas et al. (2014), Finkelstein et al. (2015), Hagen et al. (2016), and
Santos et al. (2020) found LAEs to lie above the relation, while other studies
have interpreted them as lying directly on the low-mass end of the relation
(e.g. Kusakabe et al. 2018). Figure 14 shows that the LAEs in our sample lie
largely on the SFMS, though a significant fration lie below the relation of
Sanders et al. (2018) for $M_{\star}<10^{9}\ \mathrm{M_{\odot}}$.

Figure 14: The stellar mass - star formation rate correlation for LAEs in our
sample. The trend fit by Sanders et al. (2018) for $z\sim 2.3$ star-forming
galaxies is drawn in dashed purple for comparison. LAEs in our sample largely
fall on the SFMS, though the lowest mass sources ($M_{\star}<10^{9}\
\mathrm{M_{\odot}}$) tend to fall below the relation.

In Appendix A, we explore the model-dependence of our measured galaxy
properties, since the parameters derived from SED fitting can be
systematically different using different models (see Conroy 2013). We conclude
that our results, including the median physical properties and the
correlations with $W_{\lambda}({\mathrm{Ly\alpha}})$ are not driven by our
specific choice of model.

##  5 Discussion

###  5.1 Are HETDEX LAEs Special?

The question, ‘‘What is a HETDEX LAE?’’ holds particular importance for
astronomers studying galaxy science with this survey. A vast sample of HETDEX
LAEs is upcoming, and samples of such objects selected by emission line
detection from a blind spectroscopic survey remain rare in the literature
(with the exception of the HETDEX Pilot Survey (Adams et al. 2011, Blanc et
al. 2011), which probed a smaller area to a brighter flux limit, and MUSE
surveys, which probe much smaller areas to fainter flux limits with only a
small overlap in redshift with HETDEX). Characterizing any idiosyncrasies in
the HETDEX LAE population will put these objects in context relative to the
numerous LAEs found by previous studies, and it will aid the interpretation of
future blind spectroscopic surveys for these objects in the EoR.

As described above, in our fLyα $\gtrsim 6\times 10^{-17}\ \mathrm{erg\
s^{-1}\ cm^{-2}}$ flux-limited sample (Gebhardt et al., 2021), the median
galaxy mass of $0.8^{+2.9}_{-0.5}\times 10^{9}\ \mathrm{M_{\odot}}$ lies very
close to many LAE samples selected through narrow band imaging. For example,
Gawiser et al. (2007) found a median mass of $1^{+0.6}_{-0.4}\times 10^{9}\
\mathrm{M_{\odot}}$ with a flux limit of $1.5\times 10^{-17}\mathrm{erg\
s^{-1}\ cm^{-2}}$ at $z=3.1$. Guaita et al. (2011) pushed to an even lower
median mass of $\sim 4\times 10^{8}\ \mathrm{M_{\odot}}$, roughly a factor of
two less massive than this sample’s median, with a flux limit of $2.0\times
10^{-17}\ \mathrm{erg\ s^{-1}\ cm^{-2}}$ at $z=2.1$. The MUSE HUDF went even
deeper, finding sources at $z>3$ with Ly$\alpha$ line fluxes as small as $\sim
2\times 10^{-18}\ \mathrm{erg\ s^{-1}\ cm^{-2}}$ and obtaining a median sample
mass of $\sim 2.5\times 10^{8}\ \mathrm{M_{\odot}}$. The sample of Santos et
al. (2020) was limited by medium-band line flux limits spanning $3.0-4.8\times
10^{-17}\ \mathrm{erg\ s^{-1}\ cm^{-2}}$ over $2<z<6$ (Sobral et al., 2018)
and measured a median LAE mass of $~{}2\times 10^{9}\ \mathrm{M_{\odot}}$,
consistent with this study. Of course, the mass range probed by HETDEX falls
far below samples selected using the Lyman/Lyman-alpha break (for example, the
lowest mass probed by Papovich et al. (2001) was $10^{10}\ \mathrm{M_{\odot}}$
at $2.0<z<3.5$). Thus, the HETDEX flux limit explores an LAE mass range
comparable to NB surveys, yet slightly more massive than the deepest NB and
spectroscopic surveys. At the expense of sensitivity, the HETDEX survey can
find fairly low-mass LAEs over a large continuous redshift interval, reducing
the effects of cosmic variance compared to NB observations.

As mentioned in §4.1, the LAEs in this sample do not stand out from NB samples
at similar redshift in terms of age, star formation rate, and dust extinction.
Thus, we can conclude that the HETDEX survey selects a typical LAE having
properties consistent with the general NB-selected population, but it may have
slightly higher stellar mass based on the line flux limit of the survey.

Nonetheless, our sample may stand out in its relation to the SFR-$M_{\star}$
relation shown in Figure 14. Compared to the relation measured in Sanders et
al. (2018), LAEs in the sample with $M_{\star}\lesssim 10^{9}\
\mathrm{M_{\odot}}$ appear to lie below the trend. This contrasts markedly
with the work of Hagen et al. (2016), who compiled their sample using the
HETDEX Pilot survey (Adams et al. 2011, Blanc et al. 2011) and found their
LAEs to lie above the SFMS. Interestingly, the LAEs lying below the SFMS in
Figure 14 have very high $W_{\lambda}({\mathrm{Ly\alpha}})$, which correlates
with lower $M_{\star}$ and SFR in Figure 13. We are not surprised that the
lowest mass systems in our sample have the highest values of
$W_{\lambda}({\mathrm{Ly\alpha}})$ given the negative correlation with
$M_{\star}$ and the fact that low mass objects need large
$W_{\lambda}({\mathrm{Ly\alpha}})$ to be detected by HETDEX, but their
position below the SFMS is peculiar. It could be related to the weak negative
correlation we found between $W_{\lambda}({\mathrm{Ly\alpha}})$ and sSFR, or
could simply be an artifact of our small sample size. This motivates further
study of the positions of LAEs on the SFMS with larger samples.

###  5.2 Which Properties Drive Ly$\alpha$ Emission?

While the size of the sample analyzed in this study is small, we were still
able to extract important information linking galaxy stellar-population
properties to Ly$\alpha$ emission strength. As the number of LAEs detected by
HETDEX grows in fields with rich photometric data, such as the Spitzer-HETDEX
Exploratory Large-Area Survey (SHELA) (Papovich et al., 2016), the number of
LAEs with measured galaxy properties will grow by many orders of magnitude.
This will provide a trove of useful data for explaining why some galaxies
shine brightly in Ly$\alpha$ while others do not, as well as exploring the
effects of galaxy environment on Ly$\alpha$ emission.

We found a significant, strong negative correlation between
$W_{\lambda}({\mathrm{Ly\alpha}})$ and stellar mass in our sample (see the top
left panel of Figure 13). This trend is often theoretically attributed to low
mass, star-forming galaxies having less neutral gas to resonantly scatter the
Ly$\alpha$ photon (as well as less dust) leading to a shorter total path
length to exit the galaxy without absorption by dust (see Ando et al. 2006).
In this sample, $W_{\lambda}({\mathrm{Ly\alpha}})$ also negatively correlated
(even more strongly) with SFR, and the fact that stellar mass and star
formation rate correlate strongly with each other complicates the
interpretation of this result. Weiss et al. (2021) addressed this issue by
binning their sample of [O iii]-emitting galaxies with Ly$\alpha$ line flux
measurements from HETDEX according to stellar mass and SFR. They found mass to
better predict $f_{esc}^{\mathrm{Ly\alpha}}$ at fixed SFR than SFR did at
fixed mass.

Fascinatingly, we did not find even a weak correlation between dust extinction
and $W_{\lambda}({\mathrm{Ly\alpha}})$. This seems surprising given that many
authors have noted such a correlation and that the theoretical explanation is
inarguable: resonantly scattered Ly$\alpha$ photons can get absorbed readily
in the presence of even a small amount of dust. A partial explanation for our
sample’s behavior with $A_{V}$ could be that it consists of systems exhibiting
strong Ly$\alpha$ emission, not absorption. For example, Reddy et al. (2021)
studied systems with Ly$\alpha$ in net absorption or emission and found a
strong correlation between $W_{\lambda}({\mathrm{Ly\alpha}})$ and $E(B-V)$. If
our sample contained objects with negative $W_{\lambda}({\mathrm{Ly\alpha}})$,
perhaps those objects would reveal the correlation. Nevertheless, other
studies of only emitters ($W_{\lambda}({\mathrm{Ly\alpha}})>0$) have also
noted a trend with dust extinction, such as Marchi et al. (2019), though a
close examination of their Figure 7 shows that the negative correlation is
largely driven by weak emitters with $W_{\lambda}({\mathrm{Ly\alpha}})<10\
\mathrm{\AA}$. Our small dynamic range in $W_{\lambda}({\mathrm{Ly\alpha}})$
may obfuscate a correlation with dust extinction. This interpretation may also
be complicated by the Ly$\alpha$ photon’s ability to escape the galaxy even in
the presence of large amounts of dust. Given a clumpy ISM geometry, clumps of
gas and dust can act as mirrors to Ly$\alpha$ photons, which ‘‘bounce’’ of the
surfaces of these clumps through resonant scattering by neutral gas, while
continuum photons pass through and thus experience extinction. Gronke et al.
(2016) found that simulated Ly$\alpha$ emission lines agreed well with
observations for models with clumpy ISM geometries, and Finkelstein et al.
(2009) found that clumpy-ISM models better fit the SEDs of over half their NB-
selected sample of LAEs at $z\sim 4.5$. Vargas et al. (2014) also found their
sample of 20 NB-selected LAEs at z=2.1 favored clumpy-ISM models.

Lastly, we found a moderate correlation between
$W_{\lambda}({\mathrm{Ly\alpha}})$ and galaxy mass-weighted age. The strength
of Ly$\alpha$ emission depends on both its production through recombination in
HII regions as well as its escape through channels in the ISM with low neutral
gas covering fractions, so the interplay between these processes determines
$W_{\lambda}({\mathrm{Ly\alpha}})$. As noted by Marchi et al. (2019), who
obtained a similar result, the trend with age could arise from older systems
having experienced intense star formation in their past, where stellar winds
and radiation cleared out neutral gas and dust, leaving channels for
Ly$\alpha$ escape. Through ongoing star formation or recent bursts, these
objects can still produce Ly$\alpha$ photons, and the ISM conditions favor
their escape. For the youngest galaxies, even though the most massive,
ionizing photon-producing stars are present, it is possible that a significant
amount of dust and neutral gas has yet to be swept away, hindering the escape
of Ly$\alpha$.

##  6 Predicting Lyman-alpha emission in the Epoch of Reionization

Using our knowledge of Ly$\alpha$ emission from HETDEX galaxies situated in an
ionized IGM, we can attempt to predict the intrinsic emission strength of LAEs
at $z>7$, an era where starlight from galaxies was still actively re-ionizing
the universe.

###  6.1 An LAE Sample in the Epoch of Reionization

Our sample at $1.9<z<3.5$ provides a view of Ly$\alpha$ emission unobscured by
a significant IGM neutral fraction. By creating a predictive model that
connects global galaxy properties to their intrinsic
$W_{\lambda}({\mathrm{Ly\alpha}})$ in this pristine era, we can apply it to
LAEs in the EoR to derive their expected intrinsic
$W_{\lambda}({\mathrm{Ly\alpha}})$, then attributing any deficiency of
Ly$\alpha$ emission from objects in the EoR to an increasing neutral fraction.
This does require the assumption that the production and escape of Ly$\alpha$
photons does not evolve with redshift for fixed galaxy properties, which will
require further testing. As a pilot attempt here, we took advantage of the
sample of $z>7$ LAEs that Jung et al. (2020) found in GOODS-N to test our
ability to predict Ly$\alpha$ emission from EoR galaxies.

Using a deep, spectroscopic survey conducted with Keck/MOSFIRE, Jung et al.
(2020) found 10 $>4\sigma$ Ly$\alpha$ detections at $z>7$ among 72 high-$z$
candidate galaxies. Such objects likely reside in ionized bubbles of the IGM,
allowing the Ly$\alpha$ photon to redshift away from the resonant-frequency
therefore lowering the absorption cross-section with neutral hydrogen. These
emitters thus serve as direct tests of our understanding of the galaxy
properties that modulate Ly$\alpha$ emission strength from the ISM/CGM.

Because the photometric catalog for the GOODS-N field contains the LAEs
discovered by Jung et al. (2020), we performed the same SED analysis detailed
in section 3.1 for those objects. We again masked all photometric bands
including and blueward of Ly$\alpha$ given the object’s spectroscopic
redshift. For most of the $z>7$ LAEs, this left 3 $HST$ filters as well as
both $Spitzer$/IRAC channels. We again used Bagpipes to estimate the galaxy
properties, adopting our fiducial model (delayed-$\tau$ SFH, Calzetti et al.
1994 dust law). Figure 15 shows an example fit for an object at $z=7.51$.

Figure 15: An example Bagpipes SED fit for LAE ID z7_GND_42912 at $z=7.51$
detected by Jung et al. (2020). For scaling purposes, we do not show the
upper-limits for non-detections in the $HST$ bands blueward of the Ly$\alpha$
break. From our photometric data, we constrained the stellar population
properties of ten LAEs in the EoR, allowing us to predict their intrinsic
Ly$\alpha$ emission using our HETDEX sample.

###  6.2 A Predictive Model for $W_{\lambda}({\mathrm{Ly\alpha}})$

To predict the Ly$\alpha$ equivalent widths of the $z>7$ sample, we chose
several properties that strongly impact the emergent Ly$\alpha$ emission from
galaxies: stellar mass, dust extinction, and star formation rate. As discussed
above, stellar mass may determine the amount of neutral hydrogen gas (and thus
dust) in the galaxy as well as the total path length needed to escape. In the
presence of dust, Ly$\alpha$ photons may terminate their resonant scattering
process through absorption by a dust grain following re-emission at longer
wavelengths, limiting likelihood of escape. Finally, the global star formation
rate impacts the production of UV photons that can create Ly$\alpha$ through
recombination, and feedback from star formation may impact the structure of
the ISM itself, creating ionized channels for escape.

Using the posterior distributions sampled by Bagpipes, we matched each $z>7$
emitter to LAEs in the HETDEX sample based on SED-derived properties. To do
this, we calculated the ‘‘separation" in the log mass, SFR, dust attenuation
parameter space from the EoR LAEs to each LAE in the HETDEX sample. For the
separation calculation, we divided each parameter value by the full range of
values in the sample to normalize the parameter space. For example, for log
stellar mass, an object in the HETDEX sample with log mass halfway between the
sample minimum and maximum would have a value of 0.5, so the difference
between 0.5 and the EoR LAE log stellar mass scaled the same way would become
input to the Euclidean distance formula. We then ranked the HETDEX LAEs by
separation in parameter space and constructed the prediction using the
$N=3,5,\ \mathrm{and}\ 7$ closest neighbors. We computed the posterior
$W_{\lambda}({\mathrm{Ly\alpha}})$ distribution by co-adding Gaussian
distributions with mean and standard deviation set by the
$W_{\lambda}({\mathrm{Ly\alpha}})$ measurements and error bars in our sample.
To give more importance to those LAEs that closely resembled the EoR galaxy,
we weighted each Gaussian distribution by the inverse of its squared distance
in parameter space from the EoR galaxy when co-adding to obtain the final
prediction. The predicted $W_{\lambda}({\mathrm{Ly\alpha}})$ distributions are
normalized such that the integral over all equivalent widths equals unity.

Figure 16: Probabilistic predictions of $W_{\lambda}({\mathrm{Ly\alpha}})$ for
ten LAEs at $z>7$ having emission line $S/N>4$. The distributions were
normalized by setting their integrals to unity. We chose stellar mass, SFR,
and dust extinction as predictive properties for this calculation. Gold,
magenta, and indigo lines show the probability distribution of our predictions
using N=3,5, and 7 nearest neighbors, and the gray shaded region shows the 68%
confidence interval for the equivalent width measurements from Jung et al.
(2020). Object IDs and redshifts are indicated with text for each plot. We
find good agreement between prediction and observation for the majority of
strong emitters

Figure 16 shows our predicted $W_{\lambda}({\mathrm{Ly\alpha}})$ distributions
for LAEs in the Jung et al. (2020) sample with Ly$\alpha$ $S/N>4$. We show
predictions using three different values of $N$, the number of nearest
neighbors in parameter space, to reveal any stochasticity in the prediction.
The measured Ly$\alpha$ equivalent widths from Jung et al. (2020) are
indicated by vertical dashed lines with $1\sigma$ error intervals shaded grey.
Importantly, we only expect our predictions to match the observed equivalent
widths of EoR LAEs if they exist in ionized bubbles. If the EoR LAEs instead
exist in regions of the IGM with significant neutral fractions, we expect to
over-predict the Ly$\alpha$ emission. On the other hand, an under-prediction
of the Ly$\alpha$ emission from an EoR object would imply our sample size is
too small to account for the diversity in physical properties of the LAE
population.

Figure 17: Predicted vs. observed $W_{\lambda}({\mathrm{Ly\alpha}})$ for the
LAE sample in Figure 16. We computed the first moment and square root of the
second moment of each $N=5$ distribution in Figure 16 for the predicted values
and their error bars. A right-pointing arrow indicates the predicted value of
object z7_GND_34204. A one-to-one dashed line is drawn to guide the eye, and
points above this line (the region shaded red) could be the result of IGM
absorption.

In Figure 17, we plot the predicted versus observed equivalents widths with a
one-to-one line drawn to facilitate comparison. Each object’s predicted value
and error were calculated as the first moment and square root of the second
moment of the $N=5$ curves in Figure 16, respectively. In five out of ten
cases (ID z7_GND_18626, z7_GND_44088, z7_GND_42912, z7_GND_22233, and
z7_GND_39781) the $1\sigma$ interval of our $W_{\lambda}({\mathrm{Ly\alpha}})$
predictions overlapped with the $1\sigma$ interval of the observational
measurement, indicating moderate agreement. For strong emitters (observed
$W_{\lambda}({\mathrm{Ly\alpha}})>20\ \mathrm{\AA}$), our prediction
overlapped with observation five out of eight times. Furthermore, two strong
emitters (z7_GND_42912 and z7_GND_16863), postulated by Jung et al. (2020) to
inhabit ionized bubbles, had observed equivalent widths greater than or equal
to the majority of our predicted $W_{\lambda}({\mathrm{Ly\alpha}})$
distributions, as one might expect for sources with little IGM attenuation.

It is not surprising that our model failed to predict weak Ly$\alpha$ emission
accurately. First, our model predicts Ly$\alpha$ EWs in the absence of IGM
absorption, thus an under prediction could imply significant absorption of
Ly$\alpha$ photons by neutral hydrogen in the IGM. Second, as our sample by
construction contains far more strong emitters than weak ones (see Figure 12),
this could presently bias us towards an over-prediction of Ly$\alpha$ emission
strengtdrasticallyh. We note that we under-predicted the emission from ID
z7_GND_34204 (indicated by an arrow in Figure 17), which could be attributed
to the dearth of objects in our sample with very high equivalent widths to
match with that object’s value, $\sim 280\ \mathrm{\AA}$.

ID z7_GND_42912 offers a good example of how challenging predicting Ly$\alpha$
emission can be. As $N$ increases, the peak of the predicted distribution
shifts from agreeing well with the observation to under-predicting it. It is
clear that our sample is presently too small to fully span the parameter space
in both $W_{\lambda}({\mathrm{Ly\alpha}})$ and physical properties. Future
analyses with much larger samples made possible by HETDEX should be able to
better capture the mean trends as well as variance in galaxy parameters that
determine Ly$\alpha$ emission strength.

Some of the predictions in Figure 16 bode well for constraining the expected
$W_{\lambda}({\mathrm{Ly\alpha}})$ given a suite of galaxy properties measured
from broadband SED fitting. With larger samples that suffer less from the
inherent idiosyncratic behavior of Ly$\alpha$ emission (for example, its
dependence on the observer’s line-of-sight), a rigorous, statistical
understanding of the properties that drive that emission will arise, unlocking
the potential of LAEs to probe cosmic reionization. We further note that, with
larger samples, machine learning (ML) may prove an invaluable tool in making
the nuanced connection between global galaxy properties and Ly$\alpha$
emission strength, as the problem requires a regression analysis well suited
for ML techniques.

##  7 Summary

We used SED fitting to study the properties of a sample of LAEs from the
HETDEX survey in GOODS-N to better understand the phenomenology behind
Ly$\alpha$ emission and ultimately leverage these beacons of light in the
distant Universe as probes of cosmic reionization.

To build the sample, we inspected 842 emission line detections to determine if
the line was Ly$\alpha$ or a feature from a low-redshift galaxy, such as [O
ii]. We then created a procedure to synthesize information about angular
separation from the emission line detection position, extracted emission line
flux, and $\chi^{2}$ of SED fit assuming $z_{\mathrm{Ly\alpha}}$ to identify
the continuum counterpart in our deep, mult-band $HST$ imaging in GOODS-N.
After removing detections with no counterparts, AGN contaminants, and sources
with insufficient photometric data, we analyzed a sample of 72 LAEs using SED
fitting performed by Bagpipes.

Our sample’s properties were consistent with studies of LAEs from NB imaging
surveys at similar redshifts. Our median sample mass was
$0.8^{+2.9}_{-0.5}\times 10^{9}\ \mathrm{M_{\odot}}$, and the galaxies’ SFRs
appeared to put them approximately on the star-forming main sequence, except
for at $M_{\star}<10^{9}\ \mathrm{M_{\odot}}$. Using Ly$\alpha$ emission line
flux measurements from HETDEX, we also studied correlations between
$W_{\lambda}({\mathrm{Ly\alpha}})$ and galaxy properties. We found strong
correlations between $W_{\lambda}({\mathrm{Ly\alpha}})$ and stellar mass as
well as SFR. We additionally found a moderate correlation where galaxies with
older stellar populations had larger Ly$\alpha$ equivalent widths.
Interestingly, we did not find a significant impact of dust extinction on
$W_{\lambda}({\mathrm{Ly\alpha}})$, whereas many other studies have. Overall,
this paints a picture of LAEs as low-mass systems with moderate star formation
activity wherein Ly$\alpha$ photons can escape even in the presence of dust.
Also, the LAEs detected by HETDEX do not stand out significantly in terms of
their stellar population properties from LAEs found using NB imaging with
comparable flux limits.

Finally, we used our LAE sample to try to predict the value of
$W_{\lambda}({\mathrm{Ly\alpha}})$ for ten LAEs at $z>7$ by matching the
distinct samples in the parameter space of mass, SFR, and dust extinction. Our
prediction matched the data at the $1\sigma$ level five out of ten times (5/8
for strong emitters); the three over-predictions could indicate significant
absorption by a neutral hydrogen in the IGM. With large sample sizes in the
near future and tools such as machine learning, we are optimistic about the
ability of HETDEX LAEs to unlock the potential of Ly$\alpha$ as a reliable
reionization probe.

##  8 Acknowledgements

APM and SLF acknowledge support from the National Science Foundation, through
grants AST-1908817 and AST-1614798. I.J. acknowledges support from NASA under
award number 80GSFC21M0002.

HETDEX is led by the University of Texas at Austin McDonald Observatory and
Department of Astronomy with participation from the Ludwig-Maximilians-
Universität München, Max-Planck-Institut für Extraterrestrische Physik (MPE),
Leibniz-Institut für Astrophysik Potsdam (AIP), Texas A&M University,
Pennsylvania State University, Institut für Astrophysik Göttingen, The
University of Oxford, Max-Planck-Institut für Astrophysik (MPA), The
University of Tokyo and Missouri University of Science and Technology. In
addition to Institutional support, HETDEX is funded by the National Science
Foundation (grant AST-0926815), the State of Texas, the US Air Force (AFRL
FA9451-04-2- 0355), and generous support from private individuals and
foundations.

The observations were obtained with the Hobby-Eberly Telescope (HET), which is
a joint project of the University of Texas at Austin, the Pennsylvania State
University, Ludwig-Maximilians-Universität München, and Georg-August-
Universität Göttingen. The HET is named in honor of its principal benefactors,
William P. Hobby and Robert E. Eberly.

VIRUS is a joint project of the University of Texas at Austin, Leibniz-
Institut für Astrophysik Potsdam (AIP), Texas A&M University (TAMU), Max-
Planck-Institut für Extraterrestrische Physik (MPE), Ludwig-Maximilians-
Universität München, Pennsylvania State University, Institut für Astrophysik
Göttingen, University of Oxford, and the Max-Planck-Institut für Astrophysik
(MPA).

The authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin for providing high performance computing,
visualization, and storage resources that have contributed to the research
results reported within this paper. URL: http://www.tacc.utexas.edu

The Institute for Gravitation and the Cosmos is supported by the Eberly
College of Science and the Office of the Senior Vice President for Research at
Pennsylvania State University.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543, doi: 10.1088/0067-0049/182/2/543
  * Acquaviva et al. (2011) Acquaviva, V., Gawiser, E., & Guaita, L. 2011, ApJ, 737, 47, doi: 10.1088/0004-637X/737/2/47
  * Acquaviva et al. (2012) Acquaviva, V., Vargas, C., Gawiser, E., & Guaita, L. 2012, ApJ, 751, L26, doi: 10.1088/2041-8205/751/2/L26
  * Adams et al. (2011) Adams, J. J., Blanc, G. A., Hill, G. J., et al. 2011, ApJS, 192, 5, doi: 10.1088/0067-0049/192/1/5
  * Ando et al. (2006) Ando, M., Ohta, K., Iwata, I., et al. 2006, ApJ, 645, L9, doi: 10.1086/505652
  * Ashby et al. (2015) Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2015, ApJS, 218, 33, doi: 10.1088/0067-0049/218/2/33
  * Atek et al. (2014) Atek, H., Kunth, D., Schaerer, D., et al. 2014, A&A, 561, A89, doi: 10.1051/0004-6361/201321519
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Blanc et al. (2011) Blanc, G. A., Adams, J. J., Gebhardt, K., et al. 2011, ApJ, 736, 31, doi: 10.1088/0004-637X/736/1/31
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582, doi: 10.1086/174346
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379, doi: 10.1093/mnras/sty2169
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718, doi: 10.1086/309250
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Cowie & Hu (1998) Cowie, L. L., & Hu, E. M. 1998, AJ, 115, 1319, doi: 10.1086/300309
  * Dijkstra (2014) Dijkstra, M. 2014, PASA, 31, e040, doi: 10.1017/pasa.2014.33
  * Du et al. (2018) Du, X., Shapley, A. E., Reddy, N. A., et al. 2018, ApJ, 860, 75, doi: 10.3847/1538-4357/aabfcf
  * Erb et al. (2012) Erb, D. K., Quider, A. M., Henry, A. L., & Martin, C. L. 2012, ApJ, 759, 26, doi: 10.1088/0004-637X/759/1/26
  * Farrow et al. (2021) Farrow, D. J., Sánchez, A. G., Ciardullo, R., et al. 2021, MNRAS, 507, 3187, doi: 10.1093/mnras/stab1986
  * Feltre et al. (2020) Feltre, A., Maseda, M. V., Bacon, R., et al. 2020, A&A, 641, A118, doi: 10.1051/0004-6361/202038133
  * Feroz & Skilling (2013) Feroz, F., & Skilling, J. 2013, in American Institute of Physics Conference Series, Vol. 1553, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 32nd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. U. von Toussaint, 106--113, doi: 10.1063/1.4819989
  * Finkelstein et al. (2015) Finkelstein, K. D., Finkelstein, S. L., Tilvi, V., et al. 2015, ApJ, 813, 78, doi: 10.1088/0004-637X/813/1/78
  * Finkelstein et al. (2010) Finkelstein, S. L., Papovich, C., Giavalisco, M., et al. 2010, ApJ, 719, 1250, doi: 10.1088/0004-637X/719/2/1250
  * Finkelstein et al. (2009) Finkelstein, S. L., Rhoads, J. E., Malhotra, S., & Grogin, N. 2009, ApJ, 691, 465, doi: 10.1088/0004-637X/691/1/465
  * Finkelstein et al. (2012) Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164, doi: 10.1088/0004-637X/756/2/164
  * Finkelstein et al. (2019) Finkelstein, S. L., D’Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36, doi: 10.3847/1538-4357/ab1ea8
  * Finkelstein et al. (2021) Finkelstein, S. L., Bagley, M., Song, M., et al. 2021, arXiv e-prints, arXiv:2106.13813.  https://arxiv.org/abs/2106.13813
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gawiser et al. (2007) Gawiser, E., Francke, H., Lai, K., et al. 2007, ApJ, 671, 278, doi: 10.1086/522955
  * Gebhardt et al. (2021) Gebhardt, K., Mentuch Cooper, E., Ciardullo, R., et al. 2021, arXiv e-prints, arXiv:2110.04298.  https://arxiv.org/abs/2110.04298
  * Giavalisco et al. (2004) Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93, doi: 10.1086/379232
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35, doi: 10.1088/0067-0049/197/2/35
  * Gronke et al. (2016) Gronke, M., Dijkstra, M., McCourt, M., & Oh, S. P. 2016, ApJ, 833, L26, doi: 10.3847/2041-8213/833/2/L26
  * Gronwall et al. (2007) Gronwall, C., Ciardullo, R., Hickey, T., et al. 2007, ApJ, 667, 79, doi: 10.1086/520324
  * Guaita et al. (2010) Guaita, L., Gawiser, E., Padilla, N., et al. 2010, ApJ, 714, 255, doi: 10.1088/0004-637X/714/1/255
  * Guaita et al. (2011) Guaita, L., Acquaviva, V., Padilla, N., et al. 2011, ApJ, 733, 114, doi: 10.1088/0004-637X/733/2/114
  * Hagen et al. (2016) Hagen, A., Zeimann, G. R., Behrens, C., et al. 2016, ApJ, 817, 79, doi: 10.3847/0004-637X/817/1/79
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hathi et al. (2016) Hathi, N. P., Le Fèvre, O., Ilbert, O., et al. 2016, A&A, 588, A26, doi: 10.1051/0004-6361/201526012
  * Hayes (2015) Hayes, M. 2015, PASA, 32, e027, doi: 10.1017/pasa.2015.25
  * Hayes et al. (2014) Hayes, M., Östlin, G., Duval, F., et al. 2014, ApJ, 782, 6, doi: 10.1088/0004-637X/782/1/6
  * Hill et al. (2008) Hill, G. J., Gebhardt, K., Komatsu, E., et al. 2008, Astronomical Society of the Pacific Conference Series, Vol. 399, The Hobby-Eberly Telescope Dark Energy Experiment (HETDEX): Description and Early Pilot Survey Results, ed. T. Kodama, T. Yamada, & K. Aoki, 115
  * Hill et al. (2021) Hill, G. J., Lee, H., MacQueen, P. J., et al. 2021, arXiv e-prints, arXiv:2110.03843.  https://arxiv.org/abs/2110.03843
  * Horne (1986) Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801
  * Huang et al. (2021) Huang, Y., Lee, K.-S., Shi, K., et al. 2021, arXiv e-prints, arXiv:2104.11354.  https://arxiv.org/abs/2104.11354
  * Jung et al. (2018) Jung, I., Finkelstein, S. L., Livermore, R. C., et al. 2018, ApJ, 864, 103, doi: 10.3847/1538-4357/aad686
  * Jung et al. (2020) Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2020, ApJ, 904, 144, doi: 10.3847/1538-4357/abbd44
  * Kelz et al. (2014) Kelz, A., Jahn, T., Haynes, D., et al. 2014, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, VIRUS: assembly, testing and performance of 33,000 fibres for HETDEX, 914775, doi: 10.1117/12.2056384
  * Khostovan et al. (2021) Khostovan, A. A., Malhotra, S., Rhoads, J. E., et al. 2021, MNRAS, 503, 5115, doi: 10.1093/mnras/stab778
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36, doi: 10.1088/0067-0049/197/2/36
  * Kornei et al. (2010) Kornei, K. A., Shapley, A. E., Erb, D. K., et al. 2010, ApJ, 711, 693, doi: 10.1088/0004-637X/711/2/693
  * Kusakabe et al. (2018) Kusakabe, H., Shimasaku, K., Ouchi, M., et al. 2018, PASJ, 70, 4, doi: 10.1093/pasj/psx148
  * Lang et al. (2016) Lang, D., Hogg, D. W., & Mykytyn, D. 2016, The Tractor: Probabilistic astronomical source detection and measurement.  http://ascl.net/1604.008
  * Lee et al. (2010) Lee, S.-K., Ferguson, H. C., Somerville, R. S., Wiklind, T., & Giavalisco, M. 2010, ApJ, 725, 1644, doi: 10.1088/0004-637X/725/2/1644
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Leung et al. (2017) Leung, A. S., Acquaviva, V., Gawiser, E., et al. 2017, ApJ, 843, 130, doi: 10.3847/1538-4357/aa71af
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Malhotra & Rhoads (2004) Malhotra, S., & Rhoads, J. E. 2004, ApJ, 617, L5, doi: 10.1086/427182
  * Marchi et al. (2019) Marchi, F., Pentericci, L., Guaita, L., et al. 2019, A&A, 631, A19, doi: 10.1051/0004-6361/201935495
  * Martin et al. (2015) Martin, C. L., Dijkstra, M., Henry, A., et al. 2015, ApJ, 803, 6, doi: 10.1088/0004-637X/803/1/6
  * Matthee et al. (2016) Matthee, J., Sobral, D., Oteo, I., et al. 2016, MNRAS, 458, 449, doi: 10.1093/mnras/stw322
  * Matthee et al. (2021) Matthee, J., Sobral, D., Hayes, M., et al. 2021, MNRAS, 505, 1382, doi: 10.1093/mnras/stab1304
  * Miralda-Escudé (1998) Miralda-Escudé, J. 1998, ApJ, 501, 15, doi: 10.1086/305799
  * Neufeld (1991) Neufeld, D. A. 1991, ApJ, 370, L85, doi: 10.1086/185983
  * Ono et al. (2021) Ono, Y., Itoh, R., Shibuya, T., et al. 2021, ApJ, 911, 78, doi: 10.3847/1538-4357/abea15
  * Ouchi et al. (2020) Ouchi, M., Ono, Y., & Shibuya, T. 2020, ARA&A, 58, 617, doi: 10.1146/annurev-astro-032620-021859
  * Oyarzún et al. (2017) Oyarzún, G. A., Blanc, G. A., González, V., Mateo, M., & Bailey, John I., I. 2017, ApJ, 843, 133, doi: 10.3847/1538-4357/aa7552
  * Papovich et al. (2001) Papovich, C., Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620, doi: 10.1086/322412
  * Papovich et al. (2011) Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123, doi: 10.1111/j.1365-2966.2010.17965.x
  * Papovich et al. (2016) Papovich, C., Shipley, H. V., Mehrtens, N., et al. 2016, ApJS, 224, 28, doi: 10.3847/0067-0049/224/2/28
  * Partridge & Peebles (1967) Partridge, R. B., & Peebles, P. J. E. 1967, ApJ, 147, 868, doi: 10.1086/149079
  * Pentericci et al. (2009) Pentericci, L., Grazian, A., Fontana, A., et al. 2009, A&A, 494, 553, doi: 10.1051/0004-6361:200810722
  * Pentericci et al. (2010) Pentericci, L., Grazian, A., Scarlata, C., et al. 2010, A&A, 514, A64, doi: 10.1051/0004-6361/200913425
  * Ramsey et al. (1994) Ramsey, L. W., Sebring, T. A., & Sneden, C. A. 1994, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 2199, Spectroscopic survey telescope project, ed. L. M. Stepp, 31--40, doi: 10.1117/12.176221
  * Reddy et al. (2021) Reddy, N. A., Topping, M. W., Shapley, A. E., et al. 2021, arXiv e-prints, arXiv:2108.05363.  https://arxiv.org/abs/2108.05363
  * Rhoads et al. (2000) Rhoads, J. E., Malhotra, S., Dey, A., et al. 2000, ApJ, 545, L85, doi: 10.1086/317874
  * Robertson et al. (2015) Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19, doi: 10.1088/2041-8205/802/2/L19
  * Sanders et al. (2018) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99, doi: 10.3847/1538-4357/aabcbd
  * Santos et al. (2020) Santos, S., Sobral, D., Matthee, J., et al. 2020, MNRAS, 493, 141, doi: 10.1093/mnras/staa093
  * Scarlata et al. (2009) Scarlata, C., Colbert, J., Teplitz, H. I., et al. 2009, ApJ, 704, L98, doi: 10.1088/0004-637X/704/2/L98
  * Shapley et al. (2003) Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65, doi: 10.1086/373922
  * Shimakawa et al. (2017) Shimakawa, R., Kodama, T., Shibuya, T., et al. 2017, MNRAS, 468, 1123, doi: 10.1093/mnras/stx091
  * Sobral et al. (2018) Sobral, D., Santos, S., Matthee, J., et al. 2018, MNRAS, 476, 4725, doi: 10.1093/mnras/sty378
  * Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289, doi: 10.1088/0004-637X/717/1/289
  * Trainor et al. (2019) Trainor, R. F., Strom, A. L., Steidel, C. C., et al. 2019, ApJ, 887, 85, doi: 10.3847/1538-4357/ab4993
  * Vargas et al. (2014) Vargas, C. J., Bish, H., Acquaviva, V., et al. 2014, ApJ, 783, 26, doi: 10.1088/0004-637X/783/1/26
  * Weiss et al. (2021) Weiss, L. H., Bowman, W. P., Ciardullo, R., et al. 2021, ApJ, 912, 100, doi: 10.3847/1538-4357/abedb9
  * Wold et al. (2014) Wold, I. G. B., Barger, A. J., & Cowie, L. L. 2014, ApJ, 783, 119, doi: 10.1088/0004-637X/783/2/119
  * Xue et al. (2016) Xue, Y. Q., Luo, B., Brandt, W. N., et al. 2016, ApJS, 224, 15, doi: 10.3847/0067-0049/224/2/15
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513

##  Appendix A Model-Dependence of Measured Galaxy Properties

Bayesian approaches to SED fitting, like the one implemented in Bagpipes
provide robust constraints on the parameter uncertainties and their
interdependence, but the model chosen for comparison to the data (as well as
the chosen priors) determines the accuracy of those estimates. In other words,
an inaccurate model yields inaccurate measurements of galaxy properties. Many
galaxy SED fitting studies have shown that model choices, such as the SFH,
systematically impact the measured galaxy properties (see Conroy 2013 for
review).

To test the robustness of our results to different modeling choices, we
performed an additional analysis of our entire sample using an alternate
model. We did not seek to find a more (or less) accurate model; we simply
wanted a different model to determine if the median properties or correlations
between Ly$\alpha$ emission and galaxy properties changed. To this end, we
adopted a constant SFH parametrization as well as the dust absorption model of
Charlot & Fall 2000. The constant SFH required two parameters: the time when
star formation began and the constant star formation rate. For dust
attentuation, we adopted the recipe given in Charlot & Fall (2000) by using an
absorption curve proportional to $\lambda^{-0.7}$, and a factor of three
reduction in the dust extinction normalization for stellar populations older
than $10^{7}$ years to account for the dispersal of stellar birth clouds. The
authors found this recipe to match the absorption of stellar continuum and
nebular emission for nearby starburst galaxies very well, and the differential
extinction toward young stars differs markedly from the treatment by Calzetti
et al. (1994) used in our ‘‘fiducial’’ model presented above.

Figure 18 shows the distribution of LAE properties measured using the
alternate model compared with the fiducial model. The sample median stellar
mass increased by 0.1 dex, as did the median SFR. These two changes do not
affect our results or interpretation significantly. The median dust dropped
from $A_{V}=0.30$ to 0.17, a fairly substantial change, but not unusual given
the common factors of $\sim$ a few discrepancies between different models and
SED-fitting codes (see Leja et al. 2017). Nonetheless, the correlations
between galaxy properties and $W_{\lambda}({\mathrm{Ly\alpha}})$ remained
unaffected by the model modifications, as shown in Figure 19. Stellar mass and
SFR correlated strongly and negatively with Ly$\alpha$ emission strength,
while other parameters, like dust extinction, continued to show no significant
correlations.

Figure 18: Comparison of galaxy properties as measured using our ‘‘fiducial’’
model (light blue) versus our ‘‘alternate’’ model (sea green). The 16th, 50th,
and 84th percentiles calculated using the alternate model are indicated by
vertical dashed grey lines, and their values are indicated with text. The
distributions are consistent, save for dust extinction ($\mathrm{A_{V}}$),
which has lower values by a factor of $\sim 2$ for the alternate model.
Figure 19: Comparison of correlations between
$W_{\lambda}({\mathrm{Ly\alpha}})$ and galaxy properties as measured using our
‘‘fiducial’’ model (light blue) versus our ‘‘alternate’’ model (sea green).
Binned values from the alternate model are indicated as open gold circles. The
correlations presented in §4.3 appear robust when different models are
adopted.

##  Appendix B Imaging, Emission lines, and SED fits for LAEs in this study

In this section, for all 72 LAEs in our sample, we present $HST$ imaging
cutouts in Figure 22 showing the sources and any neighbors, the HETDEX
Ly$\alpha$ emission line detections in Figure 21, and the SED fits with
Bagpipes (Carnall et al., 2018) used to measure physical properties in Figure
20.

Figure 20: All SED fits for LAEs in the sample (see Figure 9 for a description
of the plots). The $\chi^{2}$ value for each fit is also given with text.
Figure 21: All emission line detections from the HETDEX Survey for LAEs in the
sample. The observed data are indicated by the blue lines with grey error
bars. A Gaussian model fit to the data is shown in green. The x-axis is scaled
in Angstroms relative to the line center. Detection IDs and Ly$\alpha$ line
redshifts are indicatd with text. Figure 22: All $HST$ F160W ($H$-band) images
of LAEs in the sample. Each cutout shows a $5\arcsec\times 5\arcsec$ image
centered on each galaxy in our sample. The galaxy centroid is indicated with a
pink diamond, and the HETDEX detection position and PSF FWHM are indicated by
a light blue cross and dashed circle, respectively. We also include object IDs
and redshifts with text.