# A selection of H$\rm\alpha$ emitters at $z=2.1$–$2.5$ using the $K_{s}$-band
photometry of ZFOURGE

Yasunori Terao  Institute of Astronomy, The University of Tokyo, Osawa 2-21-1,
Mitaka, Tokyo 181-0015, Japan  Lee R. Spitler School of Mathematical and
Physical Sciences, Macquarie University, Sydney, NSW 2109, Australia  Research
Centre in Astronomy, Astrophysics & Astrophotonics, Macquarie University,
Sydney, NSW 2109, Australia  Australian Astronomical Optics, Faculty of
Science and Engineering, Macquarie University, Macquarie Park, NSW 2113,
Australia  Kentaro Motohara  Advanced Technology Center, National Astronomical
Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan  Department
of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1
Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Nuo Chen  Department of Astronomy,
Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan

(Accepted October 31, 2022)

###### Abstract

Large and less-biased samples of star-forming galaxies are essential to
investigate galaxy evolution. H$\rm\alpha$ emission line is one of the most
reliable tracers of star-forming galaxies because its strength is directly
related to recent star formation. However, it is observationally expensive to
construct large samples of H$\rm\alpha$ emitters by spectroscopic or narrow-
band imaging survey at high-redshifts. In this work, we demonstrate a method
to extract H$\rm\alpha$ fluxes of galaxies at $z=2.1$–$2.5$ from $K_{s}$
broad-band photometry of ZFOURGE catalog. Combined with 25–39 other filters,
we estimate the emission line fluxes by SED fitting with stellar population
models that incorporate emission-line strengths. 2005 galaxies are selected as
H$\rm\alpha$ emitters by our method and their fluxes show good agreement with
previous measurements in the literature. On the other hand, there are more
H$\rm\alpha$ luminous galaxies than previously reported. The discrepancy can
be explained by extended H$\rm\alpha$ profiles of massive galaxies and a
luminosity dependence of dust attenuation, which are not taken into account in
the previous work. We also find that there are a large number of low-mass
galaxies with much higher specific star formation rate (sSFR) than expected
from the extrapolated star formation main sequence. Such low-mass galaxies
exhibit larger ratios between H$\rm\alpha$ and UV fluxes compared to more
massive high sSFR galaxies. This result implies that a “starburst” mode may
differ among galaxies: low-mass galaxies appear to assemble their stellar mass
via short-duration bursts while more massive galaxies tend to experience
longer-duration ($>10\ \mathrm{Myr}$) bursts.

††journal: ApJ

##  1 Introduction

Since galaxies increase their masses through star formation, it is important
to investigate properties of star-forming galaxies for understanding galaxy
evolution. How actively galaxies form stars at a given epoch is typically
measured with the cosmic star formation rate density (SFRD), which is the star
formation rate (SFR) per unit volume. This rate peaked at cosmological
redshifts of $z\sim 2$–$3$ (Madau & Dickinson, 2014).

A correlation between galaxy stellar mass and SFR is also an important
statistical property of galaxy populations. This so-called star formation main
sequence (SFMS) has been observationally confirmed at least up to $z\sim 4$
(e.g., Tomczak et al., 2016). Its normalization increases with redshift, which
indicates that galaxies at higher redshift tend to have higher SFR per stellar
mass (specific star formation rate; sSFR) than in the local universe. Thus
measurements of the SFMS tell us how star formation activities in galaxies
change through cosmic time.

Scatter around the SFMS can constrain galaxy evolution models because it
likely reflects the dominant mechanism of star formation in galaxies. For
example, Tacchella et al. (2016) have suggested that observed small scatter
($\sim 0.3\ \mathrm{dex}$) can be explained by a “self-regulated” process,
where SF in galaxies is governed by a balance between supply and depletion of
gas due to inflows and star formation, respectively, rather than bursty events
such as major mergers. Other studies have also shown that the variations in
sSFR can be explained by the time evolution of gas accretion rate onto the
galaxies (Dutton et al., 2010; Dekel et al., 2013). Hence the galaxies on the
SFMS are expected to evolve along the tight correlation until they start to
quench at the high-mass end.

Sparre et al. (2017) found that simulated galaxies in the FIRE cosmological
hydrodynamical simulation (Hopkins et al., 2014) with $\log{(M/M_{\odot})}<9$
show larger scatter around the SFMS than more massive simulated galaxies. They
have attributed this finding to the increasing impact of the supernova
feedback in galaxies with shallower gravitational potentials, where outflows
caused by supernovae suppress star formation and some time later significant
inflows result in a star formation burst. Since different SFR indicators are
sensitive to star formation operating over different timescales, Sparre et al.
(2017) point out that the H$\alpha$/UV luminosity ratio is a useful metric for
capturing time since a significant change in the star formation rate of a
galaxy. Their simulated low-mass galaxies show relatively large scatter in the
H$\alpha$/UV luminosity ratios, indicating a range of star formation
timescales occurring in this galaxy population. Interestingly, recent
observations have found that galaxies both in the local universe and at $z>4$
with larger H$\alpha$/UV (Emami et al., 2019; Faisst et al., 2019). However,
these works had limited numbers as it is observationally expensive to
construct large sample of high-redshift low-mass galaxies with both the
H$\alpha$ and UV luminosities.

When measuring the cosmic SFRD and SFMS, a large survey volume is required to
reduce effects of cosmic variance and put meaningful constraints on theories
of galaxy formation. Furthermore it is important to use a relatively unbiased
selection for the sample and obtain robust SFRs. While there are various SFR
indicators available, H$\alpha$ emission line is the most direct and reliable
probe. It is produced by ionizing photons from the most massive and short-
lived stars and sensitive to short-period fluctuations of star-formation
activity. In addition, Oteo et al. (2015) have found that star-forming
galaxies traced by H$\alpha$ emission lines are sensitive to a diverse range
of galaxies at $z\sim 2$. Therefore we can obtain less-biased sample of star-
forming galaxies with robust SFR measurements by using H$\alpha$ emission line
as a tracer.

There have been various studies to probe star-forming galaxies at $z>1$ with
H$\alpha$ emission line by both spectroscopy and narrow-band imaging. For
example, Reddy et al. (2018) measured H$\alpha$ fluxes of galaxies at
$1.4<z<3.8$ and investigated their properties in combination with other
optical emission lines using spectroscopic data of the MOSDEF survey. Sobral
et al. (2014) used narrow-band imaging data of HiZELS (Geach et al., 2008;
Sobral et al., 2009) to derive H$\alpha$ luminosity functions at four redshift
bins of $z=0.40,0.87,1.47$ and $2.23$.

One limitation of spectroscopic sample is that due to their sensitivity
limits, there are very limited number of observations for low-mass
($\log{(M/M_{\odot})}\lesssim$9.0) galaxies making it challenging to make
statements about this class of galaxy (e.g., Sanders et al., 2021). An issue
with narrow-band imaging surveys is that they generally lack information about
dust attenuation. One way to overcome these limitations is to use multi-
wavelength broad-band photometry, which allows the identification of low-mass
galaxies over relatively large survey volumes. Indeed, several works have
identified H$\alpha$ emission line at $z>4$ from color excesses in broad-band
Spitzer/IRAC photometry (e.g., Smit et al., 2014, 2015; Stark et al., 2013;
Schaerer & de Barros, 2009; Faisst et al., 2016, 2019). At lower redshift,
Saito et al. (2020) and Onodera et al. (2020) have extracted H$\alpha$ fluxes
at $0.3<z<2.0$ and [OIII] fluxes at $3.0<z<3.7$, respectively, from broad-band
photometry in combination with SED fitting to multi-wavelength data. See also
Forrest et al. (2018) who extracted H$\alpha$ fluxes from composite SEDs
consisting of stacked broad-band photometry SEDs.

In this work, we estimate H$\alpha$ emission line fluxes of galaxies at $z\sim
2.1$–$2.5$ from $K_{s}$ broad-band photometry using the multi-wavelength
catalog from the ZFOURGE survey (Straatman et al., 2016). As fluxes measured
in broad-band filters contain both continuum and emission lines, it is
important to separate each contribution. We estimate continuum fluxes by SED
fitting with the multi-band photometry so that other properties of galaxies,
such as stellar mass and amount of dust extinction, can be simultaneously
estimated if the effects of emission lines are properly taken into account.

This paper is organized as follows. Section 2 describes the data and
methodology, where H$\alpha$ fluxes extracted from broad-band photometry are
compared to those from narrow-band imaging and spectroscopy. In Section 3, we
investigate properties of H$\alpha$ emitters at $z=2.1$–2.5. Especially we
discuss their H$\alpha$ luminosity function focusing on the most luminous
emitters and the low-mass end of SFMS where no large spectroscopic sample
exists. We summarize our results in Section 4

In the following we adopt the AB magnitude system and a cosmology with
$H_{0}=70\ \mathrm{km\ s^{-1}\ Mpc^{-1}},\Omega_{m}=0.3$, and
$\Omega_{\Lambda}=0.7$.

##  2 Data and method

###  2.1 ZFOURGE catalog

We use a photometric catalog of the FourStar galaxy evolution survey (ZFOURGE;
Straatman et al. 2016) for our analysis. Observations of the ZFOURGE were
conducted over 45 nights with the FourStar near-infrared camera on the
Magellan telescope (Persson et al., 2013), targeting three legacy survey
fields: COSMOS (Scoville et al., 2007), UDS (Lawrence et al., 2007), and CDFS
(Giacconi et al., 2002). The survey area is $\sim 407$ $\mathrm{arcmin^{2}}$
in total. Since a number of ground-based and space observations have been
conducted in these fields, there is significant public data so that the
catalog covers from 0.8–8 $\mathrm{\mu m}$ with 26–40 filters. The resulting
well-sampled galaxy SEDs makes it possible to accurately derive properties of
galaxies by SED fitting. The Spitzer/IRAC at 3.6 and 4.5 $\mathrm{\mu m}$
bands are key for estimate stellar masses of the galaxies by putting good
constraints on rest-frame near-infrared regimes of the SEDs at $z>1$. In
addition, the Spitzer/MIPS 24 $\mathrm{\mu m}$ band data is available for all
of the three fields, which can be used to estimate SFRs from infrared
bolometric luminosities.

A unique advantage of ZFOURGE is the very deep $K_{\mathrm{s}}$-band imaging
whose 5$\sigma$ limiting magnitude is $\geq 25.5$ across the fields, resulting
in 80% completeness for galaxies with $\log{(M/M_{\odot})}\sim 9.0$ at $z=2$
(Straatman et al., 2016).

Another unique characteristic of the ZFOURGE survey is its medium-band filters
dividing $J$ and $H$-bands into three
($J_{\mathrm{1}},J_{\mathrm{2}},J_{\mathrm{3}}$) and two ($H_{\mathrm{s}}$,
$H_{\mathrm{l}}$), respectively. These filters efficiently sample the Balmer
break (3636 Å) and 4000 Åbreak of galaxies at $1.5\lesssim z\lesssim 3.5$ and
improve photo-$z$ accuracy in this redshift range. Straatman et al. (2016)
have derived photo-$z$s from SED fitting by EAZY (Brammer et al., 2008) and
compared them with publicly available spectroscopic redshifts. They have found
median $\sigma_{z}\sim 0.013$–$0.022$ in the three fields, where
$\sigma_{z}\equiv
1.48\times|z_{\mathrm{phot}}-z_{\mathrm{spec}}|/(1+z_{\mathrm{spec}})$.

###  2.2 Flux excess in the $K_{\mathrm{s}}$-band

H$\alpha$ emission lines (rest-frame 6563 Å) of galaxies at $2.1<z<2.5$ fall
into the $K_{\mathrm{s}}$-band of the ZFOURGE, while [OIII] (rest-frame 4959 &
5007 Å) and H$\beta$ (rest-frame 4861 Å) into the two medium-bands
($H_{\mathrm{s}}$ and $H_{\mathrm{l}}$), respectively as shown in Figure 1.

Figure 1: Combinations of the emission lines and the filters at each redshift.
At $z>1$, H$\alpha$, H$\beta$, and [OIII] emission lines are observed in the
near-infrared regime. Especially, at $2.1<z<2.5$, H$\alpha$ falls into the
$K_{\mathrm{s}}$ broad-band filter of the ZFOURGE, while H$\beta$ and [OIII]
are simultaneously observed in the $H$ medium-band filters.

To derive the emission line fluxes, we define “flux excess” (in units of
$\mathrm{erg\ s^{-1}\ cm^{-2}\ }$Å-1) as a difference between an observed flux
and a continuum flux in a medium/broad-band filter as follows:

$\displaystyle f_{\textrm{excess}}=f_{\textrm{obs}}-f_{\textrm{cont}},$ (1)

where $f_{\textrm{obs}}$ and $f_{\textrm{cont}}$ are the observed flux and the
continuum flux, respectively. The continuum flux is estimated by SED fitting
taken into account of contributions of the emission lines.

###  2.3 SED fitting with emission line templates

When redshifted emission lines fall into certain filters, fluxes observed in
them are boosted by the lines. Hence the SED fitting may overestimate
continuum fluxes of the galaxies unless the effect of the emission lines is
properly taken into account. This is especially important for high-$z$
galaxies because they have been known to show larger equivalent widths
compared to those in the local universe (Stark et al., 2013; Salmon et al.,
2015; Forrest et al., 2018).

Therefore, we run an SED fitting code, Fitting and Assessment of Synthetic
Templates (FAST; Kriek et al. 2009), with templates including emission lines
to accurately estimate the continuum fluxes of the galaxies. These templates
are created by injecting emission line fluxes into a spectral templates of
Bruzual & Charlot (2003) in a manner of Salmon et al.
(2015)111https://github.com/BrettSalmon/Nebular. They have calculated the
emission-line strengths following a method of Inoue (2011), where strengths of
the UV-optical emission lines from Ly$\alpha$ to those at $\sim 1\ \mathrm{\mu
m}$ relative to H$\beta$ are calculated varying ionization parameters,
metallicities, and hydrogen number densities using CLOUDY (Ferland et al.,
1998). H$\beta$ luminosity is estimated to be:

$\displaystyle L_{\mathrm{H\beta}}=4.78\times
10^{-13}\frac{1-f_{\mathrm{esc}}}{1+0.6f_{\mathrm{esc}}}N_{\mathrm{lyc}}$ (2)

where $f_{\mathrm{esc}}$ and $N_{\mathrm{lyc}}$ represent the escape fraction
and the production rate of the Lyman continuum photons in a galaxy,
respectively. While $N_{\mathrm{lyc}}$ is calculated for each Bruzual &
Charlot (2003) model template based on each age, there is no consensus on
$f_{\mathrm{esc}}$ and its redshift dependency. Salmon et al. (2015) have
provided calculations for two extreme cases, $f_{\mathrm{esc}}=0$ and
$f_{\mathrm{esc}}=1$. We adopt $f_{\mathrm{esc}}=0$, where all the Lyman
continuum photons are absorbed and their energy is converted to the emission
lines, to be consistent with a result obtained for Lyman break galaxies at
$z\sim 3$ (Nestor et al., 2011).

Figure 2 shows several examples of the best-fit SEDs obtained with and without
the emission line templates (blue and orange, respectively). It clearly shows
that the continuum levels, which are mainly constrained by the Spitzer/IRAC
photometry, can be significantly overestimated resulting in the overestimated
stellar masses by at most $\sim 0.8\ \mathrm{dex}$.

Figure 2: Examples of the best-fit SEDs with (blue) and without (orange) the
emission line templates, where IDs are from the ZFOURGE catalog. The
Spitzer/IRAC photometry at 3.6 and 4.5 $\mathrm{\mu m}$ bands, which sample
rest-frame $\sim 2\ \mathrm{\mu m}$, can not be fitted well unless the effects
of the emission lines are taken into account. The estimated stellar masses can
differ at most $\sim 0.8\ \mathrm{dex}$. The residuals shown just below the
best-fit SEDs are the differences between the observations and the best-fit
templates including emission lines (purple crosses). On the other hand, the
flux excesses shown in the bottom panels show the differences between the
observations and the continuum fluxes estimated from the best-fit SEDs. Red
circles represent the photometry of the $K_{\mathrm{s}}$-band, in which
H$\alpha$ emission line falls. Green circles represent the H medium-band
filters of ZFOURGE, in which H$\beta$ and/or O[III] emission lines fall.

###  2.4 Robustness against model assumptions

The SED fitting described above allows us to derive not only the continuum
levels but also the emission line strengths of galaxies. However, the line
strengths in the templates strongly depend on model assumptions, since the
correct values of $f_{\mathrm{esc}}$ and $N_{\mathrm{lyc}}$ are controversial
and may vary among galaxies. In contrast, the continuum fluxes estimated by
the SED fitting are more robust against the model assumptions. Figure 3 shows
how the continuum fluxes of our sample vary when $N_{\mathrm{lyc}}$ is doubled
in the SED fitting. Each panel shows flux difference due to $N_{\mathrm{lyc}}$
as a function of continuum flux (left) and a histogram (right).

Figure 3: (Left): Changes of the estimated continuum fluxes of our sample when
assuming a doubled production rate of the Lyman continuum photons in the SED
fitting. (Right): The distribution of the continuum offsets between the two
models.

We find a typical difference of $\lesssim 10\%$ between the continuum fluxes
derived with the different assumptions. Therefore we extract the continuum
component from the SED fitting with the emission-line templates, and subtract
them from observed fluxes to obtain H$\alpha$ fluxes. In our sample, a typical
difference between the H$\alpha$ fluxes from the flux excesses versus those
derived from direct SED fitting is $\sim 40\%$.

###  2.5 Selection of H$\alpha$ emitter at $z=2.1$–$2.5$

We select our initial candidates of H$\alpha$ emitters at $z=2.1$–$2.5$ by
requiring the flux excesses in the $K_{\mathrm{s}}$-band to be
$f_{\mathrm{excess}}\geq 2\sigma_{K_{\mathrm{s}}}$, where
$\sigma_{K_{\mathrm{s}}}$ is photometric error in the band. This selection
results in 2005 H$\alpha$ emitters (699, 632, and 674 in COSMOS, UDS, and
CDFS, respectively).

H$\alpha$ fluxes (in units of $\mathrm{erg\ s^{-1}\ cm^{-2}}$) of our sample
are derived from the flux excesses as follows:

$\displaystyle f_{\mathrm{H\alpha}}=r_{\mathrm{H\alpha}}\times
f_{\mathrm{excess}}\times\Delta\lambda,$ (3)

where $r_{\mathrm{H\alpha}}$ is a ratio of H$\alpha$ to total strengths of all
emission lines in the $K\mathrm{s}$-band, $f_{\mathrm{excess}}$ is the flux
excess, and $\Delta\lambda$ is a bandwidth of the $K\mathrm{s}$ filter. The
other major emission lines expected to fall in the $K{\mathrm{s}}$-band are
[NII] and [SII]. The constant $r_{\mathrm{H\alpha}}$ is calculated from the
relative line-strengths used in Section 2.3 adopting the photo-$z$ of
individual galaxies to determine which lines fall into the $K\mathrm{s}$-band.
A median value of $r_{\mathrm{H\alpha}}$ is $\sim$0.79 over a metallicity
range from $0.20Z_{\odot}$ to $Z_{\odot}$. It increases up to
$r_{\mathrm{H\alpha}}\sim 0.9$ for an extreme case of $Z=0.02Z_{\odot}$,
however, the impact on the final $f_{\mathrm{H\alpha}}$ estimate is still
insignificant. These values are consistent with the line ratios of star-
forming galaxies at $z\sim 2$ measured by Kewley et al. (2016) from
spectroscopy.

Uncertainties of the H$\alpha$ fluxes are estimated using photometric errors
in the $K{\mathrm{s}}$-band. The H$\alpha$ flux limit is calculated for each
galaxy based on our flux excess threshold. As we will discuss in Section 3.2,
our sample is mass-limited at $\log{(M/M_{\odot})}>9.0$.

###  2.6 H$\alpha$ flux comparisons with other methods

####  2.6.1 Narrow-band imaging

The ZFOURGE catalog contains photometry of the NB209 filter, which has been
used in the New H$\alpha$ survey (Lee et al., 2012), for galaxies in COSMOS
and CDFS. The NB209 filter has been designed to have a central wavelength of
20,987 Å and a FWHM of 208 Å to detect H$\alpha$ from galaxies at $z=2.2$.

To check the consistency of our method with a conventional narrow-band emitter
selection, we derive H$\alpha$ fluxes from $(K_{\mathrm{s}}-NB209)$ color
excesses and compare them with those derived from the $K_{\mathrm{s}}$-band
excesses. Our selection procedure follows general techniques often used in
narrow-band surveys (e.g., Sobral et al., 2009; Lee et al., 2012), which
includes applying various S/N thresholds to identify high-quality flux
measurements. As a result of these quality filters, 237 galaxies are selected
as emitter candidates. We then narrow down the candidates to be at $2.1<z<2.5$
based on their photometric redshifts, which results in 73 H$\alpha$ emitters.

We find 63/73 (86%) of the narrow-band selected galaxies are also selected by
our method and the fraction does not depend on the amount of the color excess.
This suggests that our method can reproduce most of the narrow-band selected
sample regardless of their equivalent widths.

A comparison of the H$\alpha$ fluxes derived using the narrow-band color
excesses and the $K_{\mathrm{s}}$-band excesses is shown in Figure 4.

Figure 4: Comparison between the H$\alpha$ fluxes derived from the NB209 color
excesses and the flux excesses in the $K_{\mathrm{s}}$-band. Both fluxes are
in the units of $\mathrm{erg\ s^{-1}\ cm^{-2}}$. The grey line shows the one-
to-one relation while the regions shaded in grey and light-grey indicate
factor 2 and 5 intervals, respectively. The color coding of individual
galaxies corresponds to their stellar masses.

We find that 46/63 (73%) of the H$\alpha$ emitters have consistent values
within a factor of 2 in both methods. The typical scatter around the one-to-
one relation is 0.16 dex. Although some outliers from the one-to-one relation
exist, there is no clear trend with stellar mass.

####  2.6.2 Long-slit spectroscopy

Some of the ZFOURGE galaxies have been observed by a spectroscopic survey,
ZFIRE (Nanayakkara et al., 2016). We find that 93 H$\alpha$ emitters in our
sample have ZFIRE spectroscopic observations. Figure 5 shows a comparison of
the H$\alpha$ fluxes between our method and ZFIRE, where a correction for
different aperture sizes is applied.

We find that 62/93 (67%) of the H$\alpha$ emitters have consistent fluxes in
both methods within a factor of 2. The typical scatter around the one-to-one
relation is 0.21 dex. Here the values of ZFIRE are corrected for the slit loss
assuming the aperture corrections applied to individual galaxies in the
ZFOURGE, where the mean factor of the correction is $\sim 2.7$.

####  2.6.3 Integral-field spectroscopy

We also compare our fluxes with the integral field spectroscopic (IFS) data
taken by $\mathrm{KMOS^{3D}}$ survey (Wisnioski et al., 2019). In our sample,
spectra of 86 H$\alpha$ emitters have been taken by the $\mathrm{KMOS^{3D}}$
(29, 21, and 36 in COSMOS, UDS, and CDFS, respectively). Figure 6 shows a
comparison of the H$\alpha$ fluxes derived by our method and the IFS. For
53/86 (62%) of the H$\alpha$ emitters, the H$\alpha$ fluxes are consistent
within a factor of 2. The typical scatter around the one-to-one relation is
0.23 dex.

####  2.6.4 Scatter between other methods

As shown above, scatter around the one-to-one relations between our method and
the narrow-band color excess, the long-slit spectroscopy, and the integral
field spectroscopy are 0.16 dex, 0.21 dex, and 0.23 dex, respectively. To
check whether the scatter is caused by uncertainties in our method, we compare
the H$\alpha$ fluxes between the three methods already reported in the
literature.

The scatter between the ZFIRE and $\mathrm{KMOS^{3D}}$ is 0.16 dex, while the
scatter between the ZFIRE and the narrow-band (NB209) color excess is 0.22
dex. There are only four galaxies observed in both the NB209 and the
$\mathrm{KMOS^{3D}}$. We find a mean value of their flux differences is 0.23
dex, however it is difficult to extract statistical conclusion due to the
small size of the sample.

Although the numbers of galaxies are small in the above comparisons, the
scatters are similar to those found between our method and the others.
Therefore we conclude that there is no significant difference in the
uncertainties of our H$\alpha$ flux derivation compared to other methods.

In addition, we have also applied our method to $H$ medium bands where [OIII]
(4959, 5008 $\mathrm{\AA}$) and H$\beta$ emission lines of our H$\alpha$
emitters fall into. The fluxes of the [OIII] doublet have been derived from
the $H$-band excesses assuming a H$\alpha$/H$\beta$ ratio of 2.86. When
compared to spectroscopic measurements of 3D-HST (Brammer_2012; Momcheva_2016)
they show good agreement with a typical scatter of $<0.17\ \mathrm{dex}$ (N.
Chen et al., in preparation).

Figure 5: Same as Figure 4, but for the ZFIRE and our method, where the
$K_{\mathrm{s}}$-band excesses are derived using the $1^{\prime\prime}.2$
aperture photometry. Figure 6: Same as Figure 4, but for the
$\mathrm{KMOS^{3D}}$ and our method.

###  2.7 Corrections for dust extinction

When converting the H$\alpha$ fluxes to luminosities, we apply corrections for
dust extinction. According to Nordon et al. (2013), the amount of dust
extinction at rest-frame 1600 Å is derived from infrared excess (IRX) as
follows:

$\displaystyle
A_{1600}=2.5\log\left\\{{\frac{\mathrm{SFR(IR)}}{\mathrm{SFR(UV)}}+1}\right\\},$
(4)

where $\mathrm{SFR(IR)}$ and $\mathrm{SFR(UV)}$ are SFRs measured from IR
luminosity and UV continuum luminosity, respectively. The IR luminosity is
estimated by integrating the IR spectral template of Wuyts et al. (2007) over
8–1000 $\mathrm{\mu m}$ scaled to match an observed 24 $\mathrm{\mu m}$ flux,
while the UV luminosity is estimated from flux density at 1600 Å obtained by
SED fitting by EAZY (Brammer et al., 2008), assuming the flatness of spectral
slope of UV continuum between 1216–3000 Å (Bell et al., 2005).

We then convert $A_{1600}$ to $A_{\mathrm{H\alpha}}$ assuming the extinction
curve of Calzetti et al. (2000). In this conversion, we assume that an
$f$-factor, a ratio between dust extinctions of stellar continuum and nebular
continuum to be unity. Although there are uncertainties, it has been known
that the $f$-factor increases with redshift and becomes nearly unity at $z\sim
2$ (e.g., Erb et al., 2006; Koyama et al., 2015; Kashino et al., 2017).
Moreover, Faisst et al. (2019) have found that even in the local universe,
where a typical $f$-factor is $\sim 0.44$, galaxies with high H$\alpha$
equivalent widths ($\gtrsim 100$ Å) have $f$-factors close to unity. As the
redshifts of our H$\alpha$ emitters are larger than 2 and their typical
equivalent width is larger than 100 Å (as we will see later), the assumption
of $f$-factor to be unity for our sample is reasonable.

For galaxies without detections in the Spitzer/MIPS 24 $\mathrm{\mu m}$ band
(S/N $<$ 3), we use relation between visual extinction derived from SED
fitting by FAST ($A_{V,\rm{FAST}}$) and IRX extinction obtained from stacked
MIPS 24$\mu$m fluxes ($A_{V,\mathrm{IRX}}$) as outlined in the following.
First, we divide our sample into $A_{V,\rm{FAST}}$ bins regardless of the S/N
of the MIPS detections. Median of 24 $\mathrm{\mu m}$ and 1600 Å fluxes are
calculated in each bin to obtain $A_{1600}$ via Equation 4, which is then
converted to $A_{V\rm,IRX}$ assuming the extinction curve of Calzetti et al.
(2000). Figure 7 shows a correlation between$A_{V,\rm{IRX}}$ and
$A_{V,\rm{FAST}}$, where the best-fit relation is given by

$\begin{split}A_{V,\mathrm{IRX}}=0.034A_{V,\mathrm{FAST}}^{3}-0.230A_{V,\mathrm{FAST}}^{2}+1.59A_{V,\mathrm{FAST}}.\end{split}$
(5)

We use this relation to estimate visual extinctions for galaxies without IR
detection, which is then converted to $A_{\mathrm{H\alpha}}$ using the curve
of Calzetti et al. (2000) again.

We find that there is an offset between $A_{V,\mathrm{IRX}}$ and
$A_{V,\mathrm{FAST}}$. For galaxies with $\mathrm{SFR(H\alpha)}>10$, a typical
value of $A_{\mathrm{H\alpha}}$ based on the IRX is $\sim 0.21\ \mathrm{mag}$
larger than that from the SED fitting. This discrepancy is consistent with the
underlying physical tracers of dust, which differs between these methods. In
the SED fitting, attenuation is estimated from a slope of the UV continua via
an empirical relation (IRX-$\beta$ relation; e.g., Meurer et al., 1999;
Overzier et al., 2011). If a galaxy is particularlly dusty, UV emission from
star forming regions can be absorbed by dust and hence UV observations may not
reflect the true star formation rate (e.g., Qin et al., 2019). In contrast,
the IRX method indirectly captures the absorbed and re-radiated UV emission.
Therefore the IRX can better trace the total amounts of the attenuation by
dust than the SED fitting which may underestimate the attenuation in the dusty
galaxies. This is what we find in the observational data presented in Figure
7.

Figure 7: Relation between $A_{V}$ derived from the SED fitting by FAST and
the IRX. Data are median values for each bin from all galaxies without a
restriction on S/N. Orange curve shows the best-fit relation. Grey dashed line
represent the one-to-one relation as a reference.

We derive intrinsic H$\alpha$ luminosities corrected for dust extinction and
convert them to H$\alpha$ SFRs using the calibration of Kennicutt (1998) with
a correction to the Chabrier (2003) IMF:

$\displaystyle L_{\mathrm{H\alpha}}$ $\displaystyle=$ $\displaystyle
f_{\mathrm{H\alpha}}\times 4\pi d_{\mathrm{L}}^{2}\times
10^{0.4A_{\mathrm{H\alpha}}},$ (6) $\displaystyle\mathrm{SFR(H\alpha)}$
$\displaystyle=$ $\displaystyle 7.9\times 10^{-42}L_{\mathrm{H\alpha}}\times
10^{-0.24},$ (7)

where $d_{\mathrm{L}}$ is the luminosity distance corresponding to the
redshift of a galaxy.

##  3 Properties of H$\alpha$ emitters at $z=2.1$–$2.5$

###  3.1 Locations on the UVJ diagram

Figure 8 shows locations of the 2005 galaxies selected as the H$\alpha$
emitters on the restframe $UVJ$ diagram (e.g. Wuyts et al., 2007). The
restframe UVJ data is provided in the ZFOURGE catalogs.

Figure 8: Locations of the H$\alpha$ emitters on the restframe $UVJ$ diagram
(e.g. Wuyts et al., 2007). Symbol color reflects galaxy stellar mass. Red
circles indicate objects identified as AGNs by either IR, radio, or X-ray in
the catalog of Cowley et al. (2016). Black line represents the criterion of
Spitler et al. (2014) separating the star-forming and quiescent galaxies.
Identified AGNs by either IR, radio, or X-ray are indicated by the red open
circles.

It is clear that most of the H$\alpha$ emitters are classified as star-forming
galaxies when we adopt the criterion of Spitler et al. (2014) for the
classification.

Although the fraction is very small, there are 37 H$\alpha$ emitters in the
quiescent region. We note that a few low-mass galaxies with
$M_{*}<10^{9}M_{\odot}$ (shown as bluer color in Figure 8) classified as
quiescent have poor SED fitting so that their stellar masses and H$\alpha$
excesses are likely not reliable. However, most others have well-fitted SEDs
and indeed a more conservative excess selection thresholds only reduces the
H$\alpha$ emitters in the quiescent region from 37 to 22 with
3$\sigma_{K_{\mathrm{s}}}$ excess threshold. One possibility is that our
emitter selection is sensitive enough to detect weak H$\alpha$ emission lines
from galaxies in quenching phases. Indeed, the quiescent H$\alpha$ emitters
show lower SFRs than star-forming galaxies with a similar-mass. As shown later
below, they are also located below the star formation main sequence at this
redshift range.

Another possibility is contributions from unidentified AGNs. As shown in
Figure 8 there are a few known AGN in the quiescent region, which might mean
some of these have emission-line contributions from the AGN.

Finally, it is possible that these galaxies are actually dusty star-forming
galaxies who have scattered into the quiescent region due to measurement
uncertainty in the $UVJ$ data. In any case, the existence of the quiescent
H$\alpha$ emitters does not affect the results below because they only
represent 1.8% of the galaxies in H$\alpha$ luminosity range of
$\log{L_{\rm{H\alpha}}}<43.2$.

###  3.2 Deriving equivalent width

A total rest-frame equivalent width of all emission lines in the
$K_{\mathrm{s}}$-band is calculated as follows:

$\displaystyle
EW(\mathrm{all})=\frac{f_{\mathrm{excess}}}{f_{\mathrm{cont}}}\times W/(1+z),$
(8)

where $W$ is the bandwidth of the $K_{\mathrm{s}}$-band filter. The H$\alpha$
equivalent width is obtained by:

$\displaystyle EW(\mathrm{H\alpha})=r_{\mathrm{H\alpha}}\times
EW(\mathrm{all}).$ (9)

using the $r_{\mathrm{H\alpha}}$ ratio of $H\alpha$ to total strengths of all
emission lines as described in subsection 2.5.

Next, we estimate a limiting equivalent width above which the sample is
complete. As already mentioned, our equivalent width measurement is more
sensitive to strong continuum (i.e., massive) galaxies. Therefore the limiting
equivalent width is expected to show stellar mass dependence. By the
definition of the emitter selection, a threshold of equivalent width to be
selected as the emitter can be calculated as

$\displaystyle
EW_{\mathrm{lim}}(\mathrm{H\alpha})=r_{\mathrm{H\alpha}}\times\frac{2\sigma_{K_{\mathrm{s}}}}{f_{\mathrm{cont}}}\times
W/(1+z).$ (10)

As the $K_{\mathrm{s}}$ photometry errors, $\sigma_{K_{\mathrm{s}}}$, do not
vary significantly among the galaxies and $f_{\rm cont}$ correlates with
stellar masses, $EW_{\mathrm{lim}}(\mathrm{H\alpha})$ is expected to show
strong correlation with the stellar mass. We have estimated
$EW_{\mathrm{lim}}(\mathrm{H\alpha})$ for all the H$\alpha$ emitters as shown
in Figure 9. In the data we find a correlation with stellar mass of the form:

$\displaystyle\log{EW_{\mathrm{lim}}(\mathrm{H\alpha})}=-0.46\log{(M/M_{\odot})}+6.33.$
(11)

Figure 9 also shows our data generally agrees with estimates from the
literature (Sobral et al., 2014; Reddy et al., 2018). Given this agreement, it
seems our sample is at least mass complete at $\log{(M/M_{\odot})}>9.0$, which
corresponds to the 80% mass completeness limit of ZFOURGE at $z\sim 2$. A
caveat on this conclusion is that we may still have missed emitters with
relatively low equivalent widths.

Figure 9: H$\alpha$ equivalent widths as a function of stellar mass. Limiting
equivalent widths of our sample is shown by the dashed black line as well.
Since our selection is more sensitive to galaxies with strong continua, the
limiting equivalent width depends on stellar mass. The grey and red circles
represent individual galaxies and their medians in the mass bins,
respectively. The best-fit relation obtained in the range of
$\log{(M/M_{\odot})}>9$ is shown by the black line extrapolated toward the
low-mass end. The blue line and intervals show the typical H$\alpha$
equivalent width and its scatter reported by Reddy et al. (2018) for galaxies
with $\log{(M/M_{\odot})}\simeq 9.3$–$10.7$. Our method seems sensitive enough
to detect most of the H$\alpha$ emitters with $\log{(M/M_{\odot})}>9.0$, but
only extreme ones can be detected below this limit.

In Figure 9 we also fit to only galaxies with $\log{(M/M_{\odot})}>9.0$ and
find:

$\displaystyle\log{EW(\mathrm{H\alpha})}=-0.24\log{(M/M_{\odot})}+4.54.$ (12)

This fit agrees very well with the results of previous works based on
spectroscopy (Reddy et al., 2018) and narrow-band imaging (Sobral et al.,
2014).

For galaxies with $\log{(M/M_{\odot})}<9.0$, we cannot discuss their
statistical properties due to mass incompleteness from the parent ZFOURGE
sample. Nevertheless it is interesting that there are many low-mass H$\alpha$
emitters with extremely high H$\alpha$ equivalent widths ($>500$ Å). These
values are much higher than the extrapolation of the relation at
$\log{(M/M_{\odot})}>9.0$ to lower masses. Such galaxies are expected to have
high sSFRs and we will further investigate them later.

###  3.3 Attenuation by dust

Figure 10 shows a correlation between the H$\alpha$ SFR and estimated dust
attenuation, which was derived two methods as described in Section 2.7.

Figure 10: Correlation between H$\alpha$ SFRs and the amount of the dust
attenuation estimated at the H$\alpha$ emission line, $A_{\mathrm{H\alpha}}$.
The red hexagonal tiles show the $A_{\mathrm{H\alpha}}$ distribution of our
sample, where $A_{\mathrm{H\alpha}}$ values are derived by the IRX as
described in 2.7. Median values of $A_{\mathrm{H\alpha}}$ in SFR bins derived
from the IRX and the SED fitting by FAST are shown in the red circles and the
blue squares, respectively.

Clearly the galaxies with higher SFRs (i.e., more luminous in H$\alpha$) are
more obscured. This is unsurprising given that dust in the interstellar medium
(ISM) of the galaxies is thought to be mainly produced by supernovae and AGB
stars (e.g., Morgan & Edmunds, 2003). The best-fit relation between the
H$\alpha$ SFRs and the amounts of the attenuation obtained from IRX relation
in the sample is given by

$\displaystyle
A_{\mathrm{H\alpha},{\mathrm{I}RX}}=0.95\log{\mathrm{SFR(H\alpha)}}-0.51.$
(13)

A similar relation is obtained for UV+IR SFRs using only the galaxies detected
in the Spitzer/MIPS 24 $\mathrm{\mu m}$ band with S/N$>$3.

The above relation is for the $A_{\mathrm{H\alpha}}$ derived from the IRX. The
difference is consistent with Figure 7, where $A_{V,{\mathrm{FAST}}}$ shows
systematically smaller values than $A_{V,{\mathrm{IRX}}}$.

###  3.4 H$\alpha$ luminosity function

Before discussing the luminosity function itself, we evaluate completeness of
the sample with respect to H$\alpha$ luminosity by performing a simulation
using observed relations. First we generate 100,000 mock galaxies with stellar
masses following a power-law distribution with the slope of $-$1.38 according
to Tomczak et al. (2014). Then we estimate their stellar continuum fluxes in
the $K_{\mathrm{s}}$-band using a relation between the continuum fluxes and
the stellar masses of our sample, which is obtained from the SED fitting
carried out in Section 2.3. Next, H$\alpha$ luminosities of the mock galaxies
are drawn from a normal distribution with a mean and a standard deviation
corresponding to the observed star formation main sequence (SFMS) and its
scatter (see Section 3.6). These luminosities are then converted to the
excesses, where redshifts of the mock galaxies are randomly chosen from the
observed distribution between 2.1 and 2.5. Finally we obtain their total
$K_{\mathrm{s}}$-band fluxes by adding the excesses to the stellar continuum
fluxes. We then choose galaxies which have (1) $K_{\mathrm{s}}$-band flux
larger than 5$\sigma$, and (2) the excess fluxes larger than 1$\sigma$, where
$\sigma$ is a photometric error, to be selected as H$\alpha$ emitters. By
calculating fractions of the emitters above the detection limit in each
H$\alpha$ luminosity bin, we find that our sample is almost ($\sim 96$%)
complete for galaxies with $L_{\rm H\alpha}>10^{42.25}\ \mathrm{erg\ s^{-1}}$,
where the H$\alpha$ luminosities are corrected for attenuation by dust.

Number density of H$\alpha$ emitters in Hα luminosity is calculated by
dividing the number of the emitters in each luminosity bin by the relevant
comoving volume. In our case, the total comoving volume, which is defined by
the redshift range and the survey area, is $\sim 5.51\times 10^{5}\
\mathrm{Mpc^{3}}$. Then the number density in a bin is obtained as follows:

$\displaystyle\phi(\log(L_{c}))=\frac{1}{\Delta(\log{L})}\sum_{|\log{\frac{L_{i}}{L_{c}}|}<\frac{\Delta\log{L}}{2}}\frac{1}{\Delta
V},$ (14)

where $L_{c}$ is a central value of the bin. We count the number of emitters
and divide it by the total volume ($\Delta V$) and the bin width
($\Delta\log{L}=0.25$). This formula is the same as the one used in Sobral et
al. (2009). Finally we fit a Schechter function to the number densities to
determine the best-fit H$\alpha$ luminosity function, which is given as

$\begin{split}\phi(L)dL=\ln{10}\times\phi^{*}\left(\frac{L}{L^{*}}\right)^{\alpha}e^{-(L/L^{*})}\left(\frac{L}{L^{*}}\right)d\log{L},\end{split}$
(15)

where $L^{*}$, $\alpha$, and $\phi^{*}$ are a characteristic luminosity, a
faint-end slope, and a normalization, respectively.

Figure 11 shows our H$\alpha$ luminosity function with that of Sobral et al.
(2013).

Figure 11: H$\alpha$ luminosity function of our sample with the result of
Sobral et al. (2013) for a comparison, which are shown by red circles and blue
squares, respectively. Green triangles illustrate how our result changes after
removing AGNs as cataloged by Cowley et al. (2016). Each curve shows the best-
fit Schechter function with the parameters summarized in Table 1, where
1$\sigma$ fitting uncertainties of our result are shown by red dashed curves.

The best-fit parameters are listed in Table 1.

Table 1: Best-fit parameters of H$\alpha$ luminosity functions |  | $\alpha$ |  | $\log{L^{*}}$ |  | $\log{\phi^{*}}$ |
---|---|---|---|---|---|---|---
This work |  | -1.62 |  | 43.6 |  | -3.38 |
This work (w/o AGN) |  | -1.61 |  | 43.5 |  | -3.30 |
This work ($A_{\mathrm{H\alpha}}=1$) |  | -1.53 |  | 43.3 |  | -3.11 |
Sobral et al. (2013) |  | -1.59 |  | 42.9 |  | -2.78 |

It is immediately clear that our luminosity function has an excess in the
bright-end compared to that of Sobral et al. (2013), while the faint-end
slopes are consistent. We should note that the brightest two bins of our
sample completely consist of AGNs as cataloged by Cowley et al. (2016).
Nevertheless, the excess still exists even if galaxies with known AGN are
removed.

In the following we investigate inconsistency between our H$\alpha$ luminosity
function and that of Sobral et al. (2013). We consider two possible causes.
One is related to structural properties of the galaxies, while the other is
different corrections for dust extinction.

It has been proposed that galaxies form their stars in disks around central
cores and increase their sizes, which has been called inside-out growth (e.g.,
van Dokkum et al., 2010). At $0.7<z<1.5$, Nelson et al. (2016) have stacked
imaging data of the 3D-HST and carried out Sérsic profile fitting to find
radial profiles of H$\alpha$ emission lines to be more extended than stellar
components. Moreover, they have shown that the difference of the light
profiles has become larger with increasing stellar mass. Also at $z>2$, Suzuki
et al. (2019) have revealed extended H$\alpha$ profiles with AO-assisted
imaging data. These results indicate at a fixed aperture size one potentially
might miss H$\alpha$ flux at large radii.

The flux measurements contained in the ZFOURGE catalog were originally
performed with a $1.2^{\prime\prime}$ aperture on PSF-matched images and then
corrected to total fluxes using growth curves from the $K_{\mathrm{s}}$-band
(Straatman et al., 2016). In the HiZELS, Sobral et al. (2009) used a fixed
$2^{\prime\prime}$ aperture for all galaxies so they might have missed some
fractions of the H$\alpha$ flux at radii $>2^{\prime\prime}$. To estimate how
much their fixed aperture can affect the luminosity function, we perform a
simulation as follows. First we assume that the stellar continuum profile of a
galaxy follows a Sérsic profile (Graham & Driver, 2005), that is

$\displaystyle
I(R)=I_{\mathrm{e}}\left\\{-b_{n}\left[\left(\frac{R}{R_{\mathrm{e}}}-1\right)^{1/n}\right]\right\\},$
(16)

where $n$ is a Sérsic index, $R_{\mathrm{e}}$ a half-light radius, and
$I_{\mathrm{e}}$ a surface intensity at that radius, $I(R_{\mathrm{e}})$. Then
$b_{n}$ is determined to satisfy the following equation:

$\displaystyle\Gamma(2n)=2\gamma(2n,b_{n}),$ (17)

where $\Gamma(n)$ and $\gamma(n,x)$ are complete and incomplete gamma
functions, respectively.

We then create H$\alpha$ profiles having the Sérsic indices measured by Nelson
et al. (2016), which are average values derived from the stacked images of the
3D-HST. To reproduce the observation, we convolve the profiles with a point
spread function (PSF) corresponding to the seeing. According to Trujillo et
al. (2001), the seeing can be modeled as a moffat profile:

$\displaystyle\mathrm{PSF}(r)=\frac{\beta-1}{\pi^{2}\alpha}\left[1+\left(\frac{r}{\alpha}\right)^{2}\right]^{-\beta},$
(18)

where $\beta=2.5$. A FWHM of the seeing is given using $\alpha$ and $\beta$ by

$\displaystyle\mathrm{FWHM}=2\alpha\sqrt{2^{1/\beta}-1}.$ (19)

In the case of the HiZELS, FWHM of the seeing is $\sim 1^{\prime\prime}$
(Geach et al., 2008). We create artificial images with the above models and
perform photometry with varying aperture sizes. As a result, we obtain
cumulative fractions of H$\alpha$ fluxes as a function of radius shown in
Figure 12.

Figure 12: Cumulative fractions of H$\alpha$ fluxes as a function of radius in
our simulation. They are obtained assuming the Sérsic profiles with the
structural parameters of Nelson et al. (2016) and the seeing of the HiZELS
observation. The color coding reflects stellar masses of the galaxies, while
the grey vertical dashed line shows the physical scale corresponding to the
$2^{\prime\prime}$ aperture used in Sobral et al. (2013) at a distance of
$z=2.23$.

A radius corresponding to the $2^{\prime\prime}$ aperture at $z=2.23$ is $\sim
8\ \mathrm{kpc}$ and represented by the grey dashed vertical line in Figure
12. We can clearly see that at least $\sim 20$% of the H$\alpha$ flux are
missed by their aperture. This appears to be especially the case for the most
massive galaxies, which are expected to reside in the bright-end of the
luminosity function. We therefore conclude that it is possible that the fixed
aperture used by Sobral et al. (2013) might have missed $\sim 40$% of the
H$\alpha$ flux. We should note, however, that the measurements of radial
profiles of the H$\alpha$ fluxes in Nelson et al. (2016) are truncated at a
radii of $\sim 10\ \mathrm{kpc}$, and we have extrapolated their profiles up
to 40 kpc for our simulation.

The effects of the missing fluxes can be also clearly seen in Figure 13, where
we show a luminosity function derived from ZFOURGE data without the aperture
corrections (magenta triangles). The luminosity function is suppressed as a
whole due to missing H$\alpha$ flux.

A potential caveat of this analysis is the assumption we make in about the
H$\alpha$ sizes of our $z=$2.1–2.5 sample. It is well-established that
galaxies evolve in size at least in terms of their continuum properties (e.g.,
van der Wel et al., 2014). As we are using sizes from a lower redshift from
Nelson et al. (2016), there is a chance that at $z\sim 2$ these average sizes
are not appropriate, which will have implications for the analysis of missing
H$\alpha$ flux. Indeed, Wilman et al. (2020) have found H$\alpha$ seems to
follow continuum sizes, which means our H$\alpha$ sizes we adopt for this
analysis are likely overestimated. This would act to reduce the amount of
missing H$\alpha$ flux in the fixed size aperture. However, from van der Wel
et al. (2014), we estimate the size evolution for star-forming galaxies over
$z=1.5$ to $z=2$ only changes by $\sim 10\%$, which means more than 20% of
H$\alpha$ flux could have been missed even at $z>2$. For comprehensive
understanding on the missing flux in H$\alpha$, more observations are
necessary to reveal structural properties of galaxies at $z>2$.

Figure 13: H$\alpha$ luminosity function derived from the ZFOURGE data without
the aperture corrections assuming $A_{\mathrm{H\alpha}}=1$ (orange triangles).
The difference between this work (grey circles) and Sobral et al. (2013) (blue
squares) is almost completely gone.

Another difference between this work and Sobral et al. (2013) is the dust
extinction correction. As shown in Figure 10 there is a positive correlation
between H$\alpha$ luminosity and the estimated attenuation by dust. Sobral et
al. (2013) adopted $A_{\mathrm{H\alpha}}=1\ \mathrm{mag}$ for all the
galaxies, which is not supported by our data. They therefore have likely
underestimated intrinsic luminosities of the brightest H$\alpha$ emitters and
overestimated the luminosities for the faintest galaxies. To demonstrate how
the adoption of a fixed attenuation impacts a luminosity function, we show a
luminosity function derived from ZFOURGE data without the aperture corrections
and assuming a fixed H$\alpha$ attenuation of 1 mag with orange points in
Figure 13. This modified luminosity function agrees with that of HiZELS very
well.

From the above results, we conclude that the bright-end discrepancy between
our luminosity function and that derived in Sobral et al. (2013) is caused by
intrinsically extended H$\alpha$ profiles and the larger
$A_{\mathrm{H\alpha}}$ in bright H$\alpha$ emitters. This implies there are
larger numbers of intrinsically bright H$\alpha$ emitters than previously
reported. The main caveat on this conclusion is that number of such objects is
small, therefore the impact of cosmic variance and AGNs might complicate our
interpretation. Larger surveys and spectroscopic follow-ups are required to
better understand the population in the bright-end.

###  3.5 Cosmic star formation rate density

By integrating the luminosity function, we can derive a cosmic SFRD at
$2.1<z<2.5$. Here we remove a contribution from identified AGNs, which is
$\sim 13$% of the total H$\alpha$ luminosity. In addition, we perform the
integration with two intervals because the faint-end slope of the luminosity
function is uncertain. One is the range of $L_{\mathrm{H\alpha}}=0$–$10^{45}\
\mathrm{erg\ s^{-1}}$, while the other is of
$L_{\mathrm{H\alpha}}=10^{41.6}$–$10^{45}\ \mathrm{erg\ s^{-1}}$. These
intervals are same as those used in Sobral et al. (2013), where
$L_{\mathrm{H\alpha}}=10^{41.6}\ \mathrm{erg\ s^{-1}}$ corresponds to
$0.01L^{*}$.

Given that the H$\alpha$ emission line is a direct SFR tracer and that
H$\alpha$ emitters cover diverse populations of star-forming galaxies (Oteo et
al., 2015), we expect a consistent H$\alpha$ SFRD with previous studies based
on UV and/or IR observations. Figure 14 shows our SFRD with those previously
reported as summarized in Madau & Dickinson (2014). We find that our estimate
agrees well with previous work, suggesting that our H$\alpha$ emitter
selection reflects the underlying star-forming galaxy population and that our
correction for the dust extinction properly works to recover total H$\alpha$
SFRs. However, we cannot rule out the possibility that galaxies strongly
obscured by the dust are missing from our sample because our selection
requires galaxies to be detected in the $K_{\mathrm{s}}$-band.

Figure 14: Our cosmic star formation rate density measurement compared with
previous studies. Red and magenta circles represent our measurement with
integration range of $L_{\rm H\alpha}=0-10^{45}$ (erg s-1) and
$10^{41.6}-10^{45}$ (erg s-1), respectively. Blue squares and green triangles
represent the results of Sobral et al. (2013), with same integration range as
ours. AGN contributions have been removed from both values. The grey circles
and crosses show the results of previous studies summarized in Madau &
Dickinson (2014) for a comparison, which are measured from UV and/or IR
luminosities.

The SFRD agreement is interesting given we find a larger number of bright
H$\alpha$ emitters than previous work. This reflects the relatively small
contribution of the brightest galaxies to the total star formation activity
due to their small numbers.

###  3.6 Star formation main sequence

We also investigate the relation between stellar mass and star formation rate
of our H$\alpha$ emitter sample. In the following analysis, all the AGNs
identified by either X-ray, IR, or radio (38/2005) are removed while the
quiescent H$\alpha$ emitters are included.

Figure 15 shows a stellar masses–H$\alpha$ SFRs plot of our sample with SFMS
of previous studies (Speagle et al., 2014; Whitaker et al., 2014; Behroozi et
al., 2013).

Figure 15: H$\alpha$ star formation main sequence (SFMS) of our sample. Grey
dots represent individual galaxies while black circles show median values in
each stellar mass bin. By fitting in the range of $9<\log{(M/M_{\odot})}<11$,
we obtain the best-fit SFMS shown by the black solid line. Green triangles are
the typical values among many previous studies using observations at various
wavelengths (Behroozi et al., 2013). Black dashed and dash-dotted lines
indicate the best-fit SFMSs of Speagle et al. (2014) and Whitaker et al.
(2014), respectively. Red crosses represent the quiescent H$\alpha$ emitters
with well-constrained SEDs classified by the $UVJ$ diagram. Dotted vertical
line shows the stellar mass completeness limit of our sample while dotted
horizontal line shows the SFR limit of our sample. This limit corresponding to
the flux excess threshold to be selected as a H$\alpha$ emitter.

The black line represents our best-fit SFMS, which is given by

$\displaystyle\log{\mathrm{SFR(H\alpha)}}=0.66\log{(M/M_{\odot})}-5.32.$ (20)

In the fitting, we only used the galaxies in the range of
$9<\log{(M/M_{\odot})}<11$ because the incompleteness becomes significant
below this range. At $\log{(M/M_{\odot})}>11$ the number of the galaxies is
small and a large fraction of them are classified as quiescent. Our result is
consistent with the previous studies. Especially, a good agreement with
Whitaker et al. (2014), where SFRs were derived from UV+IR luminosities,
indicates that our correction for the dust extinction successfully recovers
total H$\alpha$ SFRs. The small offset in the most massive bins, though within
the scatter, is probably due to different selections of the star-forming
galaxies. For example, selection by observed FIR fluxes may bias a sample
toward dusty starburst systems rather than normal galaxies on the SFMS
(Rodighiero et al., 2011; Lee et al., 2013). In contrast, H$\alpha$ selection
is thought to reflect a more diverse galaxy population (Oteo et al., 2015).

###  3.7 Low-mass galaxies with high sSFR

At $\log{(M/M_{\odot})}>9$, we find that the scatter around the SFMS is $\sim
0.3\ \mathrm{dex}$, which is defined as a standard deviation of SFR offsets of
individual galaxies from the SFMS. This value agrees well with previous
observations (e.g. Shivaei et al., 2015) and support a “self-regulated”
evolution as outlined in Tacchella et al. (2016).

At $\log{(M/M_{\odot})}<9$, however, there are many low-mass galaxies with
much higher SFRs than predicted from an extrapolation of the SFMS. We should
note, of course, that only strong H$\alpha$ emitters can be detected in the
low-mass end due to the sensitivity limit of the $K_{\mathrm{s}}$-band of
ZFOURGE. Nevertheless it is still interesting to investigate them further
because the high sSFRs may suggest their bursty star formation histories
(SFHs). Several works have also favor a bursty SFH interpretation for low-mass
galaxies both observationally and theoretically (e.g. Broussard et al., 2019;
Faisst et al., 2019; Emami et al., 2019; Sparre et al., 2017).

Before investigating the SFHs, we first evaluate the impact of photometric
errors on the sSFRs of the low-mass galaxies. It might be possible that the
apparently high sSFRs are only caused by increased photometric errors in their
H$\alpha$ fluxes. To explore this possibility, we perform a simulation as
follows. First we generate 100,000 mock galaxies with stellar masses of
$\log{(M/M_{\odot})}=7.0$–$10$ following the stellar mass function of Tomczak
et al. (2014). Their SFRs are calculated by Equation 20, where we assume that
our best-fit linear relation to the SFMS is held even at
$\log{(M/M_{\odot})}<9$. Then we add fluctuations to the simulated SFRs based
on the photometric errors in the $K_{\mathrm{s}}$-band (grey circles in Figure
16).

Figure 16: SFMS plot same as Figure 15 but individual observed galaxies are
colored according to luminosity ratios between H$\alpha$ and UV (see Section
3.8). Simulated observations are shown as grey dots with the scatter in
certain stellar mass bins shown as orange bars. The difference between the
simulations and observations at $\log{(M/M_{\odot})}<9$ suggest that high
photometric scatter is not the likely causing the high SFRs at low masses.

Next we derive probability distributions from the intrinsic scatter of the
simulated observation in different stellar mass bins (orange bars in Figure
16). The larger observed scatter ($\gtrsim 0.5\ \mathrm{dex}$) compared to the
simulated scatter from photometric errors at $\log{(M/M_{\odot})}<9.0$ likely
means the high SFR values are not due to photometric errors. Therefore we
conclude that the low-mass galaxies with the high sSFRs show starburst-like
sSFRs.

###  3.8 H$\alpha$/UV ratio

If a galaxy has undergone an instantaneous starburst within the past 100
$\mathrm{Myr}$, its luminosity ratio between H$\alpha$ and UV is expected to
be different from a galaxy that is forming stars at a constant rate (Sparre et
al., 2017; Faisst et al., 2019). This is because an instantaneous burst is
reflected effectively immediately in the H$\alpha$ luminosity compared to the
observed UV. As a consequence, when a short starburst occurs, the H$\alpha$/UV
ratio is boosted (relative to a constant SFH) during first 10 $\mathrm{Myr}$.
Subsequently the ratio is returns to a lower value 10–100 $\mathrm{Myr}$ after
the burst. This means that if the high sSFRs of the low-mass galaxies have
been caused by recent starbursts, they should also have high H$\alpha$/UV
ratios.

Figure 17 shows a comparison of the H$\alpha$ luminosities and the UV
luminosities of our H$\alpha$ emitters, where both have been corrected for the
dust extinction.

Figure 17: Comparison between the H$\alpha$ and UV luminosities, which are
both corrected for the dust attenuation. The color coding reflects stellar
masses of the galaxies. The grey line shows the expected relation from a
constant SFH without starbursts within $100\ \mathrm{Myr}$ (Faisst et al.,
2019). The shaded region indicates a $\pm 0.3\ \mathrm{dex}$ region about the
relation.

Grey line represents the expected relation when assuming a constant SFH
(Faisst et al., 2019) and shaded area shows a $\pm 0.3$ dex region about the
relation. Most of the galaxies with $\log{(M/M_{\odot})}<9.0$ have very high
H$\alpha$/UV ratios, which can not be explained by a constant SFH. Very few
massive galaxies show such high ratios. Figure 16 makes it clear that the low-
mass galaxies that are above the extrapolated SFMS correspond to the
population with high H$\alpha$/UV ratios. These results suggest that the low-
mass H$\alpha$ emitters with high sSFRs found here likely recently experienced
short period starbursts.

###  3.9 SFHs of the starburst galaxies

Here we bring the results together and explore the possibility that there are
two starbust modes: one for low-mass galaxies and another for high-mass
galaxies. While in the previous section we mostly focused on low-mass
starburst, there are also high-mass galaxies with high sSFRs. From Figure 16
it is clear their H$\alpha$/UV ratios are lower than those observed for the
low mass galaxy starbursts. This may imply that starburst galaxies with
different stellar masses have different SFHs. In other words, the mechanisms
which triggered the starbursts may depend on the stellar masses of the galaxy.

To investigate this idea more quantitatively, we define the starburst galaxies
as those with SFRs larger than expected from the SFMS by 0.3 dex and measure
the median H$\alpha$/UV ratio in each stellar mass bin, as shown in Figure 18.

Figure 18: H$\alpha$/UV ratios of the starburst galaxies as a function of
stellar mass. The blue circles represent median values in each stellar mass
bin. The grey horizontal line shows the value expected from a constant SFH,
while the dashed vertical line indicates the stellar mass completeness limit.

It is clear that H$\alpha$/UV decreases with increasing stellar masses and
approaches a value expected from a constant SFH. Before interpreting this
trend, it is important to understand if any selection affects the high
H$\alpha$/UV ratios of the starburst galaxies, as H$\alpha$/UV ratios can be
overestimated if UV-faint galaxies at fixed H$\alpha$ luminosities and/or
H$\alpha$-bright galaxies at fixed UV luminosities are preferentially selected
somehow. The former is unlikely because the optical data are much deeper than
the IR, so we do not explore this possibility here. The latter, a preferential
selection of H$\alpha$-bright galaxies at fixed UV luminosities is plausible,
so we explore it further in Figure 19. The Figure shows the H$\alpha$/UV
ratios of the starburst galaxies in various stellar mass bins. As the
starburst galaxies are, by definition, selected to have high H$\alpha$
luminosities (i.e. SFRs), the luminosity criteria used for starburst selection
are always higher than the limiting luminosity for HAE selection. Given this
finding, we conclude that any selection bias inherent to the sample should
have no significant impact on the H$\alpha$/UV ratios.

Figure 19: Same as Figure 17, but only the starburst galaxies are plotted
separately according to their stellar masses. Dashed black horizontal lines
represent thresholds of the H$\alpha$ SFR to be classified as the starburst
galaxies in our definition while dashed-dotted horizontal lines is our
H$\alpha$ emitter selection. As the luminosity criteria used for starburst
selection are always higher than the limiting luminosity for HAE selection, we
conclude the ratios in Figure 18 are not likely due to a selection bias.

Therefore the impact of the bias toward luminous H$\alpha$ emitters on the
H$\alpha$/UV ratios seems insignificant. We therefore conclude that the trend
in Figure 18 is real.

The results explored here seem to imply there are two types of starburst
modes. As mentioned above, a high H$\alpha$/UV ratio implies a bursty SFH,
suggesting that a starburst has occurred within past $\sim 10\ \mathrm{Myr}$.
Therefore the trend shown in Figure 18 indicates that there are at least two
different types of starbursts. One is a burst with a short duration, which is
dominant in low-mass galaxies, causing high H$\alpha$/UV ratios. The other is
a burst with a longer duration, which is dominant in high-mass galaxies,
resulting in constant SFHs lasting for $>10\ \mathrm{Myr}$. Although it is
beyond the scope of this work to reveal specific mechanisms that govern
starburst activities in each stellar mass bin, several phenomena have been
proposed as triggers such as supernova feedback, gas inflows, and major
mergers.

##  4 Summary

We have identified 2005 H$\alpha$ emitters at $2.1<z<2.5$ from the ZFOURGE
catalog. We take advantage of the deep ZFOURGE $K_{\mathrm{s}}$ photometry and
wide wavelength coverage and derive H$\alpha$ flux estimates from observed
$K_{\mathrm{s}}$ photometry. We specifically measure the $K_{\mathrm{s}}$ flux
excess which is the difference between observed fluxes and a stellar continuum
flux estimated by SED fitting. When performing the SED fitting, contributions
of emission lines were included in the SED templates.

We have corrected H$\alpha$ luminosity for dust extinction based on the IRX
method. For galaxies without IR detection, we used a relation between $A_{V}$
derived from the IRX and the SED fitting. As a result, we have found that the
attenuation of H$\alpha$ depends on both H$\alpha$ luminosity and UV+IR SFR,
ranging from 0 to 2 mag. The correction for dust extinction is more important
for galaxies with higher SFR.

The H$\alpha$ luminosity function shows an excess in the luminous end compared
to the result of Sobral et al. (2013). The excess is mostly explained by
missing H$\alpha$ flux in Sobral et al. (2013) due to their use of a fixed
photometric aperture and their not accounting for a SFR dependence on the
strength of the dust extinction. Specifically, our simulation suggests that
the $1.2^{\prime\prime}$ aperture used in Sobral et al. (2013) would have
missed $\sim 40$% of the H$\alpha$ fluxes for the luminous-end galaxies.
Moreover, they have underestimated intrinsic H$\alpha$ luminosities of the
most luminous galaxies by at most 1 mag as they assumed
$A_{\mathrm{H\alpha}}=1$ mag for all galaxies, which is likely an incorrect
assumption given what we find in Figure 10. Therefore we conclude there are
more H$\alpha$ luminous galaxies than previously reported. The impact of this
finding on the SFRD is not important – our SFRD is consistent with other works
because SFR density contribution from the most H$\alpha$-luminous galaxies is
small. Larger surveys and spectroscopic follow-ups will allow us to further
investigate this rare H$\alpha$-bright population.

The SFMS from our sample agrees well with previous results above the stellar
mass completeness limit, $\log{(M/M_{\odot})}>9$. Interestingly, we find there
are many H$\alpha$ emitters with large sSFR below this mass limit. We show
this population cannot be caused by scatter due to photometric errors. We have
measured the H$\alpha$/UV luminosity ratio of this population, which is a good
indicator of starburst timescales. We find that the low-mass galaxies with
large sSFR tend to have large ratios, suggesting that they experienced a
shorter-timescale starburst recently. In contrast, we find evidence that
higher mass starburst galaxies have lower ratios and hence likely experienced
a starbust mode for much longer timescales.

In this work, we have verified a method to derive H$\alpha$ fluxes from
excesses in broad-band photometry, which enables us to efficiently construct
large less-biased SFG sample. This method also can be applied to other
emission lines and/or other redshifts. For example, the $J$ and $H$ medium-
band filters of ZFOURGE can detect [OII] and H$\rm\beta$+[OIII] of galaxies in
the same redshift range as H$\rm\alpha$ emitters in this work, respectively.

## Acknowledgements

We thank the anonymous reviewer for helpful comments that improved the paper.
We thank Ivo Labbé for discussions related to this work. This work was
supported by JSPS Overseas Challenge Program for Young Researchers and
MEXT/JSPS KAKENHI Grant Number 15H02062, 18H03717, and 20H00171. LS
acknowledges support from an Australian Research Council Discovery Project
grant (DP190102448).

## References

  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Bell et al. (2005) Bell, E. F., Papovich, C., Wolf, C., et al. 2005, ApJ, 625, 23, doi: 10.1086/429552
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Broussard et al. (2019) Broussard, A., Gawiser, E., Iyer, K., et al. 2019, ApJ, 873, 74, doi: 10.3847/1538-4357/ab04ad
  * 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
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Cowley et al. (2016) Cowley, M. J., Spitler, L. R., Tran, K.-V. H., et al. 2016, MNRAS, 457, 629, doi: 10.1093/mnras/stv2992
  * Dekel et al. (2013) Dekel, A., Zolotov, A., Tweed, D., et al. 2013, MNRAS, 435, 999, doi: 10.1093/mnras/stt1338
  * Dutton et al. (2010) Dutton, A. A., Van Den Bosch, F. C., & Dekel, A. 2010, MNRAS, no, doi: 10.1111/j.1365-2966.2010.16620.x
  * Emami et al. (2019) Emami, N., Siana, B., Weisz, D. R., et al. 2019, ApJ, 881, 71, doi: 10.3847/1538-4357/ab211a
  * Erb et al. (2006) Erb, D. K., Steidel, C. C., Shapley, A. E., et al. 2006, ApJ, 647, 128, doi: 10.1086/505341
  * Faisst et al. (2019) Faisst, A. L., Capak, P. L., Emami, N., Tacchella, S., & Larson, K. L. 2019, arXiv:1909.03076 [astro-ph].  http://arxiv.org/abs/1909.03076
  * Faisst et al. (2016) Faisst, A. L., Capak, P., Hsieh, B. C., et al. 2016, ApJ, 821, 122, doi: 10.3847/0004-637X/821/2/122
  * Ferland et al. (1998) Ferland, G., Korista, K., Verner, D., et al. 1998, PASP, 110, 761, doi: 10.1086/316190
  * Forrest et al. (2018) Forrest, B., Tran, K.-V. H., Broussard, A., et al. 2018, ApJ, 863, 131, doi: 10.3847/1538-4357/aad232
  * Geach et al. (2008) Geach, J. E., Smail, I., Best, P. N., et al. 2008, MNRAS, 388, 1473, doi: 10.1111/j.1365-2966.2008.13481.x
  * Giacconi et al. (2002) Giacconi, R., Zirm, A., Wang, J., et al. 2002, ApJS, 139, 369, doi: 10.1086/338927
  * Graham & Driver (2005) Graham, A. W., & Driver, S. P. 2005, PASA, 22, 118, doi: 10.1071/AS05001
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581, doi: 10.1093/mnras/stu1738
  * Inoue (2011) Inoue, A. K. 2011, MNRAS, 415, 2920, doi: 10.1111/j.1365-2966.2011.18906.x
  * Kashino et al. (2017) Kashino, D., Silverman, J. D., Sanders, D., et al. 2017, ApJ, 835, 88, doi: 10.3847/1538-4357/835/1/88
  * Kennicutt (1998) Kennicutt, J. 1998, ARA&A, 36, 189, doi: 10.1146/annurev.astro.36.1.189
  * Kewley et al. (2016) Kewley, L. J., Yuan, T., Nanayakkara, T., et al. 2016, ApJ, 819, 100, doi: 10.3847/0004-637x/819/2/100
  * Koyama et al. (2015) Koyama, Y., Kodama, T., Hayashi, M., et al. 2015, MNRAS, 453, 879, doi: 10.1093/mnras/stv1599
  * Kriek et al. (2009) Kriek, M., van Dokkum, P. G., Labbe, I., et al. 2009, ApJ, 700, 221, doi: 10.1088/0004-637X/700/1/221
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599, doi: 10.1111/j.1365-2966.2007.12040.x
  * Lee et al. (2012) Lee, J. C., Ly, C., Spitler, L., et al. 2012, PASP, 124, 782, doi: 10.1086/666528
  * Lee et al. (2013) Lee, N., Sanders, D. B., Casey, C. M., et al. 2013, ApJ, 778, 131, doi: 10.1088/0004-637X/778/2/131
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Meurer et al. (1999) Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64, doi: 10.1086/307523
  * Morgan & Edmunds (2003) Morgan, H. L., & Edmunds, M. G. 2003, MNRAS, 343, 427, doi: 10.1046/j.1365-8711.2003.06681.x
  * Nanayakkara et al. (2016) Nanayakkara, T., Glazebrook, K., Kacprzak, G. G., et al. 2016, ApJ, 828, 21, doi: 10.3847/0004-637X/828/1/21
  * Nelson et al. (2016) Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016, ApJ, 828, 27, doi: 10.3847/0004-637X/828/1/27
  * Nestor et al. (2011) Nestor, D. B., Shapley, A. E., Steidel, C. C., & Siana, B. 2011, ApJ, 736, 18, doi: 10.1088/0004-637X/736/1/18
  * Nordon et al. (2013) Nordon, R., Lutz, D., Saintonge, A., et al. 2013, ApJ, 762, 125, doi: 10.1088/0004-637X/762/2/125
  * Onodera et al. (2020) Onodera, M., Shimakawa, R., Suzuki, T. L., et al. 2020, ApJ, 904, 180, doi: 10.3847/1538-4357/abc174
  * Oteo et al. (2015) Oteo, I., Sobral, D., Ivison, R. J., et al. 2015, MNRAS, 452, 2018, doi: 10.1093/mnras/stv1284
  * Overzier et al. (2011) Overzier, R. A., Heckman, T. M., Wang, J., et al. 2011, ApJ, 726, L7, doi: 10.1088/2041-8205/726/1/L7
  * Persson et al. (2013) Persson, S. E., Murphy, D. C., Smee, S., et al. 2013, PASP, 125, 654, doi: 10.1086/671164
  * Qin et al. (2019) Qin, J., Zheng, X. Z., Wuyts, S., Pan, Z., & Ren, J. 2019, arXiv:1909.13505 [astro-ph].  http://arxiv.org/abs/1909.13505
  * Reddy et al. (2018) Reddy, N. A., Shapley, A. E., Sanders, R. L., et al. 2018, ApJ, 869, 92, doi: 10.3847/1538-4357/aaed1e
  * Rodighiero et al. (2011) Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40, doi: 10.1088/2041-8205/739/2/L40
  * Saito et al. (2020) Saito, S., de la Torre, S., Ilbert, O., et al. 2020, MNRAS, 494, 199–217, doi: 10.1093/mnras/staa727
  * Salmon et al. (2015) Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183, doi: 10.1088/0004-637X/799/2/183
  * Sanders et al. (2021) Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19, doi: 10.3847/1538-4357/abf4c1
  * Schaerer & de Barros (2009) Schaerer, D., & de Barros, S. 2009, A&A, 502, 423, doi: 10.1051/0004-6361/200911781
  * Scoville et al. (2007) Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1, doi: 10.1086/516585
  * Shivaei et al. (2015) Shivaei, I., Reddy, N. A., Shapley, A. E., et al. 2015, ApJ, 815, 98, doi: 10.1088/0004-637x/815/2/98
  * Smit et al. (2014) Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58, doi: 10.1088/0004-637X/784/1/58
  * Smit et al. (2015) Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122, doi: 10.1088/0004-637X/801/2/122
  * Sobral et al. (2014) Sobral, D., Best, P. N., Smail, I., et al. 2014, MNRAS, 437, 3516, doi: 10.1093/mnras/stt2159
  * Sobral et al. (2013) Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128, doi: 10.1093/mnras/sts096
  * Sobral et al. (2009) Sobral, D., Best, P. N., Geach, J. E., et al. 2009, MNRAS, 398, 75, doi: 10.1111/j.1365-2966.2009.15129.x
  * Sparre et al. (2017) Sparre, M., Hayward, C. C., Feldmann, R., et al. 2017, MNRAS, 466, 88, doi: 10.1093/mnras/stw3011
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15, doi: 10.1088/0067-0049/214/2/15
  * Spitler et al. (2014) Spitler, L. R., Straatman, C. M. S., Labbé, I., et al. 2014, ApJ, 787, L36, doi: 10.1088/2041-8205/787/2/L36
  * Stark et al. (2013) Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129, doi: 10.1088/0004-637X/763/2/129
  * Straatman et al. (2016) Straatman, C. M. S., Spitler, L. R., Quadri, R. F., et al. 2016, ApJ, 830, 51, doi: 10.3847/0004-637X/830/1/51
  * Suzuki et al. (2019) Suzuki, T. L., Minowa, Y., Koyama, Y., et al. 2019, PASJ, 71, 69, doi: 10.1093/pasj/psz047
  * Tacchella et al. (2016) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 457, 2790, doi: 10.1093/mnras/stw131
  * Tomczak et al. (2014) Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85, doi: 10.1088/0004-637X/783/2/85
  * Tomczak et al. (2016) —. 2016, ApJ, 817, 118, doi: 10.3847/0004-637X/817/2/118
  * Trujillo et al. (2001) Trujillo, I., Aguerri, J., Cepa, J., & Gutiérrez, C. 2001, MNRAS, 328, 977, doi: 10.1046/j.1365-8711.2001.04937.x
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28, doi: 10.1088/0004-637x/788/1/28
  * van Dokkum et al. (2010) van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018, doi: 10.1088/0004-637X/709/2/1018
  * Whitaker et al. (2014) Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104, doi: 10.1088/0004-637X/795/2/104
  * Wilman et al. (2020) Wilman, D. J., Fossati, M., Mendel, J. T., et al. 2020, ApJ, 892, 1, doi: 10.3847/1538-4357/ab7914
  * Wisnioski et al. (2019) Wisnioski, E., Schreiber, N. M. F., Fossati, M., et al. 2019, arXiv:1909.11096 [astro-ph].  http://arxiv.org/abs/1909.11096
  * Wuyts et al. (2007) Wuyts, S., Labbe, I., Franx, M., et al. 2007, ApJ, 655, 51, doi: 10.1086/509708